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I .  INTRODUCTION 


Elastic  and  inelastic  scattering  of  ion  beams  from  solid  surfaces 
has  become  a  powerful  analysis  tool  used  to  obtain  a  variety  of  infor¬ 
mation  about  such  surfaces  and  thin  films.  Film  layer  thicknesses, 
concentrations  of  heavy  atoms  in  light  substrates  and  in  many  cases 
layer  stoichiometries  are  obtainable  from  Rutherford  backscatte^ing 
(RBS)  and  nuclear  reaction  (NR)  spectra.  Chu,  Mayer  and  Nicolet  have 
written  a  comprehensive  book  on  the  subject  in  which  many  formulae  are 
derived  relating  basic  data  such  as  energy  losses  and  cross  sections  to 
layer  thicknesses  and  concentrations.  For  the  most  part,  their  dis¬ 
cussions  apply  to  situations  where  geometries  are  simple,  layers  are 
well  defined  and  of  uniform  concentrations,  and  cross  sections  are 
Rutherford  or  constant. 


In  this  report,  a  computer  program  which  we  have  developed  to 
model  highly  complex  surfaces  is  discussed.  A  surface  which  may  have 
many  intermixed  elements  concentrated  nonuniformly  in  depth,  gradual 
interfaces  between  layers,  small  amounts  of  contaminants  covering  the 
outer  surface  or  whose  analysis  depth  is  such  that  the  probing  cross 
section  varies  strongly,  cannot  be  described  by  simple  analytical  ex¬ 
pressions.  In  such  complicated  surfaces,  all  surface  constituents  con¬ 
tribute  to  the  energy  losses  and  the  experimentally  measured  energy 
distributions  are  due  to  a  complex  folding  of  yields  from  all  constitu¬ 
ents.  In  addition,  effects  such  as  concentration  gradients,  energy  loss 
straggling,  and  surface  roughness  can  combine  with  the  above  factors 
to  complicate  the  analysis  even  further. 

The  code  generates  the  expected  scattered  particle  spectra  result¬ 
ing  from  ion  bombardment  of  the  surface  in  such  a  way  that  a  direct 
comparison  to  experimentally  measured  distributions  can  be  made.  The 
only  variables  in  the  calculation  are  the  concentration  profiles  of  all 
expected  constituents  of  the  surface.  For  this  reason,  it  is  necessary 
to  have  a  prior  knowledge  of  what  elements  make  up  the  surface.  Since 
only  spectral  shapes,  not  absolute  magnitudes,  are  required  for  the  fit, 
absolute  quantities  related  to  geometries  and  beam  charge  integrations 
need  not  be  measured.  However,  the  results  are  absolute  in  that  the 
relative  atomic  percentages  of  the  target  constituents  as  a  function  of 
the  depth  in  the  surface  are  obtained.  The  basic  data  which  are  used 
in  this  calculation  are  cross-sections  and  energy  losses  of  the  pro¬ 
jectiles  in  the  target  constituents  and  the  ionization  potentials  of 
all  elements. 


A  preliminary  version  of  this  code  was  described  in  reference  2. 

In  that  report,  it  was  shown  that  the  results  obtained  with  shape  fitting 


1.  W-K  Chu,  J.W.  Mayer  and  M.A.  Nicolet,  Backscattering  Spectrometry, 
Academic  Press,  New  York,  1978. 

2.  A.  Niiler,  J.E.  Youngblood,  S.E.  Caldwell  and  T.J.  Rock,  "An 
Accelerator  Technique  for  the  Study  of  Ballistic  Surfaces",  BRL 
Report  No.  1815,  August  1975.  (AD  #A016899) 
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were  in  excellent  agreement  with  absolute  measurements  depending  on 
precise  geometry  and  beam  charge  integration*  Many  significant  changes 
have  since  been  incorporated  into  the  calculation.  They  are  of  the 
type  which  make  the  program  more  general,  provide  up-to-date  input  data, 
allow  for  a  more  accurate  determination  cf  items  such  as  energy  loss 
straggling,  or  are  of  a  housekeeping  nature.  The  original  version 
was  operational  on  an  SEL-86  computer  but  the  current  version  runs  on 
the  ARRADCOM/BRL  CDC-7600  machine. 

In  section  II  of  this  report,  a  description  of  the  calculation  is 
given  and  the  function  of  each  of  the  subroutines  of  the  program  is 
described.  Sample  calculations  of  three  different  surfaces  are  shown  in 
section  III.  The  appendices  give  the  program  listing  as  well  as  a 
description  of  the  input  cards. 


II.  PROGRAM  DESCRIPTION 


This  program  calculates  the  yield  as  a  function  of  energy  resulting 
from  a  beam  of  hydrogen  or  helium  ions  backscattering  or  interacting  with 
the  constituents  of  a  solid  surface  or  thin  film.  A  typical  backscatter- 
ing  geometry  is  assumed  where  the  incoming  beam  is  defined  by  a  set  of 
apertures,  the  target  can  be  set  at  any  angle  with  respect  to  the  beam 
and  an  energy  dispersive  particle  detector  is  set  at  any  angle  behind  an 
aperture.  The  target  is  divided  into  N  layers  (where  N  <  100)  and  each 
of  the  layers  can  be  further  subdivided  into  slabs  if  additional  calcu- 
lational  resolution  is  desired.  Each  layer  may  contain  up  to  five  dif¬ 
ferent  elements  (this  number  is  easily  increased) ;  the  concentration  of 
each  element  is  given  in  terms  of  relative  atomic  percent.  Each  layer 
is  also  assigned  a  thickness.  The  combination  of  layer  thicknesses  and 
atomic  concentrations  for  the  whole  surface  define  the  concentration 
profiles  of  the  elements  in  the  surface.  The  objective  of  the  calcula¬ 
tion  is  to  arrive  at  a  unique  set  of  concentration  profiles  which  yields 
the  scattered  particle  energy  distribution  which  best  fits  the  experi¬ 
mentally  measured  distribution. 

The  calculation  proceeds  by  allowing  the  incident  beam  to  penetrate 
into  the  first  layer  of  the  surface.  Figure  1  shows  the  beam,  surface, 
detector  geometry  schematically.  The  energy  in  the  layer  (position  x) 
is  given  by 


(1) 


where  x^  =  x/cos<j)  and  the  energy  loss  ir.  the  layer  is  given  by  the  Bragg 
additivity  rule 
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The  f j 1 s  are  the  relative  atomic  percentages  of  the  constituents  and 

the  dE|/dx!s  are  the  corresponding  specific  energy  losses.  The  pro¬ 
jectile  then  interacts  with  one  of  the  layer  constituents  and  either 
elastically  scatters  backward  at  angle  0  or  produces  a  reaction  particle 
which  is  emitted  at  0.  The  probability  cf  this  interaction  taking  place 
is  given  by  the  cross  section  cr(E^(x),0).  The  outgoing  particle  energy 

is  given  by 


Eq(x)  =  K  E.(x)  (3) 

where  K  is  the  usual  kinematic  factor.  As  the  outgoing  particle  now 
leaves  the  layer,  it  loses  more  energy  so  that  its  final  energy  becomes 


Ef(x)  =  Eq(x) 


2  dE  (x;  ) 
o 

- 7 - 

dx 


dx; 


(4) 


where  x^  =  x/cos  (tt-0-cJ))  . 

Gaussian  energy  broadening  is  applied  at  two  points  in  this  cal¬ 
culation.  The  E^(x)  in  equation  1  is  broadened  for  beam  energy  reso¬ 
lution  and  energy  loss  straggling  due  to  the  incident  energy  loss.  The 
resulting  incident  energy  distribution  Gaussian  can  be  further  divided 
into  any  number  of  energies,  each  centered  at  equal  areas  along  the 
distribution.  This  procedure  has  the  effect  of  treating  the  E_.  (x) 

distribution  as  if  it  consisted  of  many  different  energies,  E,  around 
E^ (x) ,  each  having  a  scattering  probability  at  x  determined  by  the  cross 

section,  cr(E).  Each  component  of  E^(x)  is  now  scattered  and  followed 

out  of  the  sample  according  to  equations  3  and  4  and  the  final  energy, 
E^(x)j  is  broadened  for  energy  loss  straggling  due  to  outgoing  energy  loss 

and  detector  resolution.  The  E^(x)  Gaussians  are  now  distributed  into 

a  final  spectrum,  each  having  been  appropriately  weighted  by  cross  sec¬ 
tions,  layer  thicknesses  and  constituent  concentrations.  Similar  final 
energy  distributions  are  obtained  from  all  N  layers  of  the  surface  while 
each  layer  is  assumed  to  interact  with  tne  same  number  of  incident  parti¬ 
cles,  and  the  energies,  energy  losses  and  energy  spreads  are  tracked. 

When  the  full,  final  distribution  due  to  one  surface  constituent  element 
is  completed,  the  same  procedure  is  followed  for  the  other  elements  and 
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their  contributions  added  to  the  final  energy  distribution  at  the 
proper  energies  with  the  proper  weights.  Additional  calculations 
can  be  performed  if  the  effects  of  two-dimensional  geometries  are 
important.  This  is  usually  the  case  if  the  beam  spot  or  detector 
solid  angle  is  unusually  large. 

The  program  proceeds  to  fitting  the  calculated  distribution  to 
an  experimentally  measured  one.  The  two  are  fit  by  least  squares  methods, 
and  the  concentration  profiles  are  adjusted  and  the  full  calculation  run 
until  a  satisfactory  fit  is  achieved.  The  primary  criteria  for  accept¬ 
able  fits  are  general  agreement  between  the  calculated  and  experimental 
shapes  along  with  the  lowest  possible  chi  squared  for  the  full  calcu¬ 
lation.  The  chi  squared  is  given  by 

?  1  N  7 

x  =  ITT  l  CA-ycC±)  -  Ye(i))  /Ye(i)  (5) 

i 

where  K  is  the  number  of  points  fit,  I  is  the  number  of  degrees  of  free¬ 
dom,  is  the  calculated  distribution,  is  the  experimental  distri¬ 

bution  and  A  is  the  least  squares  proportionality  constant  between  the 
calculated  and  experimental  distributions.  Typically,  I  is  given  by 

I  =  N  (Nc  -  1)  (6) 

where  N  is  the  number  of  layers  used  in  the  calculation  and  is  the 

number  of  constituents  per  layer.  Since  usually  it  is  possible  to 
obtain  an  acceptable  fit  within  five  to  ten  iterations  of  the  profiles 
by  hand,  no  automatic  search  routine  has  been  incorporated  in  this  pro¬ 
gram.  Computer  costs  would  be  driven  up  significantly  by  such  auto-fit 
routines . 


The  primary  functions  of  each  of  the  program  elements  are: 

(1)  Program  PROFILE.  This  is  the  control  program.  Inputs  read  by  this 
section  are  the  title,  the  reaction  parameters,  geometrical  factors  and 
surface  layer  configurations.  It  sets  up  the  arrays  of  reaction  param¬ 
eters  which  determine  the  kinematic  factors  for  each  projectile  -  target 
constituent  combination,  determines  whether  the  geometry  is  valid, 
follows  the  projectile  energy  losses  both  incoming  and  outgoing,  and 
calls  the  various  routines  needed  for  some  inputs,  array  setups  and 
outputs . 


(2)  Subroutine  SETIQN.  This  routine  is  used  to  calculate  and  set  up  the 
array  of  effective  ionization  potential  Ue£f)  for  all  layers.  I  ^  is 


given  by 


ln  ^eff-1 


y  z.f.  ln  (I.) 
v  i  i  ^  iJ 


l  z.f. 

V  11 


(7) 
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where  z.,  f.  and  I.  are  the  atomic  number,  percentage  concentration 
and  ionization  potential  respectively  of  zhe  i  target  constituent. 

3 

The  reference  ionization  potentials,  1^,  are  obtained  from  J.E.  Turner  . 

(3)  Subroutine  SETTAB.  This  routine  is  used  to  input  concentration 
information,  set  up  tables  of  concentrations  of  the  various  elements  in 
all  layers  and,  when  desired,  calculate  and  output  energy  loss  reference 
data  as  well  as  composite  energy  losses  in  all  layers. 


(4)  Subroutine  WIDTH.  This  routine  calculates  the  energy  loss  straggling 
from  the  Bohr  model  using  Lindhard  and  Scharff's  correction  factor  for 
low  and  medium  energy  projectiles  (see  Reference  1).  The  straggling 
width  is  given  by 


0  W  •  dE 
2  m 

a  =  - : 


w 

2-  In  — 

eff 


(8) 


where  k  =  1  for  high  energy  projectiles  and 
1  4E 

k  =  y  —  for  intermediate  and  low  energy  projectiles.  The 
eff 

numerical  expression  defining  the  line  between  high  and  intermediate 
energies  as  well  as  definitions  of  the  ocher  terms  of  equation  6  are 
given  in  section  2.6  of  reference  1.  A  sense  switch  option  permits  the 
straggling  width  to  be  set  to  zero. 


(5)  Subroutines  PICKE,  FINDE,  IGAUSS.  These  routines  are  used  to  cal¬ 
culate  integrals  of  Gaussian  distributions. 


(6)  Subroutine  DIST.  This  routine  builds  the  final  energy  distribution 
by  summing  all  outgoing  particle  gaussian  distributions  weighted  by  the 
appropriate  reaction  cross  section,  target  constituent  concentration 

in  the  appropriate  layer  and  that  layer  thickness.  This  distribution 
is  built  only  between  some  preset  minimum  and  maximum  energies. 

(7)  Subroutine  KIN.  This  routine  calculates  the  non-relativistic  two- 
body  kinematics  between  an  incoming  and  a  particular  target  particle. 

It  also  calls  the  cross  section  routine  and  converts  the  cross  section 
obtained  into  the  laboratory  system. 


3.  J.E.  Turner,  "Studies  in  Penetration  of  Charged  Particles  in 
Matter”,  NAS-NRC  Publication  11,  Nuclear  Science  Series  Report 
No.  39,  (1962). 
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(8)  Subroutine  SIG.  All  reaction  cross  sections  for  a  given  reaction 
energy  and  angle  are  determined  in  this  routine.  If  the  reaction 
Q-value  is  zero,  it  automatically  calculates  the  Rutherford  cross  section. 
However,  in  the  case  of  proton  elastic  scattering  from  oxygen,  the 
Rutherford  cross  section  is  modified  by  a  factor  obtained  from  Chow, 

4 

et  al.  .  Otherwise,  it  must  use  stored  arrays  and  interpolates  for  the 
cross  section  of  interest.  Currently,  only  the  12C(d,p  ),  14N(d,p  ) 
and  0(d,p^)  C.M.  cross  sections  at  160°  lab  are  stored.  See  references 
5,  6  and  7  for  the  C,N  and  0  cross  sections  respectively. 

(9)  Subroutine  SHAPE.  This  routine  calculates  the  mean  energy  and 
standard  deviation  of  the  energy  as  the  incoming  or  outgoing  particles 
traverse  all  but  the  final  layer  along  their  trajectories. 


(Id)  Subroutine  ELOSS.  This  routine  calculates  the  energy  loss  of  the 
Hydrogen  or  Helium  isotopes  through  any  given  target.  It  uses  the 
analytical  expression 


where 


S  =  l  f; 


*1  + 


S|  =  Al (i)  •  EAZ(l 1 


(9) 


(10) 


and 


S 


i 

H 


A3(i) 

E 


ln(l 


.  Mil) 

E 


+  A5(i)  •  E)  . 


(ID 


The  f.’s  are  the  constituent  concentrations,  E  is  the  particle  energy 
and  the  A' s  are  coefficients  empirically  determined  by  Anderson  and 

4.  H.C.  Chow,  G.M.  Griffiths  and  T.H.  Hall,  Gan.  J.  Phys.,  53  (1975) 
1672. 

5.  M.  Huez,  L.  Quaglia  and  G.  Weber,  "Fonction  D'Excitation  de  la 

12  13 

Reaction  C(d,pQ)  C  F.ntre  400  et  1350  keV  -  Distributions  Angu- 
laires",  Nucl.  Instr.  Meth.  105  (1972)  197. 

6.  A.  Niiler  and  R.  Birkmire,  "The  ^N(d,p_f^N  Cross  Section,  0.40- 

1.45  MeV",  4th  International  Ion  Beam  Analysis  Conference,  Aarhus, 
Denmark,  June  25-29,  1979. 

7.  N.  Longequeque ,  H.  Beaumevieil le,  E.  Ligeon,  J.P.  Longequeque, 

and  M.  Sandon,  "Etude  Des  Reactions  ^0(d,a  ),  ^Ofdjp  )  et 
16  o  o 

0(d,pp  de  300  keV  a  1  MeV  (Resultats  Experimentaux) ,  Le 

Journal  de  Physique,  26_  (1965)  367. 
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Ziegler^* ^  for  the  target  constituent.  Equation  9  is  the  Bragg 
additivity  rule  for  energy  losses  of  composite  materials. 

(11)  Subroutine  OUTPUT.  This  routine  inputs  the  experimental  spectrum 
along  with  the  calibration  constants  associated  with  it.  It  condenses 
the  calculated  spectrum,  arranges  it  to  match  the  experimental  channel 
energies,  calculates  the  least  squares  fit  proportionality  constant  and 
the  chi  squared  between  the  two  spectra,  and  outputs  a  listing  of  the 
two  normalized  spectra. 

(12)  Subroutine  PL0T1.  This  routine  is  used  to  output  both  experimental 
and  calculated  spectra  on  a  line  printer  plot. 

(13)  Subroutine  SWITCH.  This  routine  is  used  to  input  sense-switches 
that  control  various  calculations  and  outputs. 

(14)  Subroutine  CQNCTAB.  This  routine  calculates  and  outputs  the  abso- 

2 

lute  concentration  (in  at/cm  )  of  the  various  target  constituents  as  a 
function  of  the  target  layer. 


III.  SAMPLE  CALCULATIONS 

In  this  section,  the  results  obtained  using  PROFILE  in  the  analysis 
of  three  different  surfaces  are  discussed.  The  first  is  a  case  where  a 
well-polished,  flat  steel  slab  (FeO-7)  was  heated  in  air  to  oxidize  the 
outer  layers.  The  second  shows  results  on  a  nozzle  shaped  steel  sur¬ 
face  (Nozzle  10)  which  was  exposed  to  burning  propellant  gases  at  high 
temperatures  (T^  ^  2400°k)  and  pressures  (170  MPa) .  These  first  two 

samples  were  prepared,  bombarded  with  0.665  MeV  deuteron  beams  from 
the  ARRADCOM/BRL  Cockcroft-Waltcn  accelerator  and  analyzed  at  BRL.  The 
third  sample,  an  anodized  aluminum  slab,  was  prepared  and  bombarded  with 

1  MeV  protons  and  2  MeV  alphas  by  Simons  et  al . ^  of  Naval  Surface 
Weapons  Center,  White  Oak,  MD.  The  proton  and  alpha  spectra  were  sent 
to  BRL  for  analysis. 

Figure  2  shows  the  concentration  profiles  for  C,N,0  and  Fe  which 
were  required  to  produce  the  fits  to  the  RBS  distribution  and  NR  peaks 
shown  in  figure  3.  The  carbon  and  oxygen  profiles  are  surface  peaked 


8.  H.H.  Anderson  and  J.F.  Ziegler,  Hydrogen  Stopping  Powers  and  Ranges 
in  All  Elements,  Pergamon  Press,  New  York,  1977. 

9.  J.F.  Ziegler,  Helium  Stopping  Powers  and  Ranges  in  All  Elements, 
Pergamon  Press,  New  York,  1977. 

10.  D.  Simons,  private  communication. 
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CONCENTRATION  PROFILES 
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Figure  2.  Concentration  profiles  for  heated  iron  (FeO-7)  slab. 
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Figure  3.  Fits  to  the  Rutherford  Backscattering  and  Nuclear  Reaction  peaks  from  sample 
FeO-7  using  profiles  shown  in  figure  2. 


but  they  also  extend  into  the  sample  the  full  0.6  mg/cirT  limit  of  the 
analysis  capability  with  the  .665  MeV  deuterons.  Upon  closer 
examination,  it  can  be  seen  that  the  carbon  and  oxygen  together  make  up 

2 

all  but  a  few  percent  of  the  outermost  10  yg/cm  layer.  The  history  of 
this  sample  helps  explain  this  profile.  The  oxygen  was  introduced  upon 
heating  but  the  surface  carbon  resulted  from  vacuum  system  carbon  being 
cracked  onto  the  steel  during  beam  bombardment.  The  fact  that  there 

2 

is  carbon  at  the  4%  level  all  the  way  to  0.6  mg/cm  and  oxygen  at  the 

2 

1%  level  to  0.4  mg/cm  is  somewhat  surprising  and  is  most  likely  due  to 
alloying  impurity  or  recoil  implantation.  This  sample  has  been  used  as 
a  test  sample  for  beam  and  electronics  set  ups  for  years.  The  original 
heating  of  the  sample  is  the  most  probable  cause  of  the  relatively  small 

2 

amount  of  nitrogen  that  is  found  down  to  150  yg/cm  .  The  iron  profile 
reflects  the  fact  that  the  total  atomic  concentrations  at  all  layers 
must  equal  100%  and  that  the  full  RBS  distribution  represents  elastic 

2 

deuteron  yield  from  the  outer  .45  mg/cm  . 

Figures  4  and  5  show  the  concentration  profiles  and  fits  respectively 
to  the  steel  nozzle  exposed  to  the  high  temperature  gases  of  burning 
propellant.  As  can  be  seen,  all  the  profiles  (C,N  and  0)  are  consider¬ 
ably  different  than  what  was  found  on  simple  heating  of  steel.  Although 

the  concentrations  for  C  and  0  at  the  surface  are  still  the  highest,  the 

surface  peaks  are  not  extremely  sharp  and  the  fall-off  to  the  interior 
is  much  more  gradual.  In  addition,  there  is  a  much  larger  concentration 
of  nitrogen.  It  appears  that  the  exposure  has  produced  layers  of  oxides, 

carbides  and  nitrides  in  the  steel  surface  extending  as  far  as  . 5  to 

2 

.6  mg/cm  into  the  surface.  The  consequences  of  such  layer  formation  to 
the  integrity  of  the  surface  upon  succeeding  exposures  is  currently  being 
investigated . 

The  final  sample  is  the  anodized  aluminum  slab.  Figure  6  shows  the 
concentration  profiles  for  H,C,0  and  A1  that  were  needed  to  obtain  the 
fits  to  the  proton  and  alpha  backscattering  spectra  shown  in  figures  7 
and  8.  The  profiles  for  three  different  calculations  are  shown  in  figure 
6 . 


The  first  calculation  shown  by  the  solid  lines  in  figures  6,  7  and 

2 

8,  gave  the  best  fit.  The  whole  5  yg/cm  outer  layer  consisted  of  car- 

2 

bon,  extending  at  lower  levels  to  almost  the  100  yg/cm  depth.  Hydrogen 
was  found  at  a  constant  4%  level  throughout  the  oxide  layer.  Below  the 
carbon  surface  layer,  the  oxygen/ aluminum  ratio  is  58/38  which  is  con¬ 
sistent  with  A190^  stoichiometry  of  60%-0,  40%-Al  along  with  some  of 

the  oxygen  being  taken  up  by  H^O.  In  order  to  show  the  effect  on  the 

fit  of  a  relatively  minor  change  in  the  profiles,  a  calculation  was 
done  with  the  profiles  indicated  by  the  dot-dashed  curves.  These  pro- 
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CONCENTRATION  PROFILES 
STEEL  EXPOSED  TO  HIGH  TEMPERATURE  GASES 
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Figure  4.  Concentration  profiles  for  steel  nozzle  (Noz  #10)  exposed  to  high  temperature 
propellant  gases. 
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Figure  6.  Concentration  profiles  for  anodized  aluminum  sample.  See  the  text  for 
explanation  of  the  different  curves. 
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Figure  7.  Fit  to  the  proton  RBS  spectrum  from  the  anodized  aluminum  sample  using 
profiles  shown  in  figure  6.  See  text  for  explanation  of  the  two  curve 


CM 


<313 IA  3AliV13a 


to 

0 

txO 

> 

c 

P 

*H 1 

P 

o 

P 

o 

0 

£ 

rH 

p 

& 

6 

0 

oj 

x: 

cn 

p 

g 

m 

p 

o 

C 

•  H 

c 

6 

o 

P 

•H 

i— ( 

P 

Oj 

Cj 

C 

T3 

Oj 

0 

1— I 

N 

PH 

•H 

X 

T3 

0 

O 

c 

P 

oj 

o 

m 

0 

JZ 

p 

p 

X 

0 

e 

p 

o 

p 

0 

P 

0 

CO 

P 

, 

P 

vD 

o 

0 

0 

Ph 

P 

t/) 

P 

co 

•H 

CQ 

m 

DS 

c 

Oj 

•  H 

X 

Ph 

C 

rH 

3: 

Ctj 

o 

0 

P 

0 

O 

rH 

P 

•H 

P 

P 

O 

•H 

p 

tL, 

PH 

00 

0 

U 

P 

W5 

•H 

Uh 


files  differ  from  the  solid  lines  only  in  the  oxide-metal  interface  and 
produce  a  discernable  change  in  the  fit  only  to  the  alpha  spectrum  and 
only  in  the  region  below  0.45  MeV.  The  combined  chi-squared  increased 
by  40%  over  the  best  fit  value.  The  final  calculation  was  done  assum¬ 
ing  no  hydrogen  at  all  in  the  oxide  layer,  shown  by  the  dashed  curves 
in  figures  6,  7  and  8.  The  alpha  spectrum  was  fit  with  a  curve  over¬ 
lapping  the  solid  line  in  all  areas,  but  this  set  of  profiles  gave  a  fit 
to  the  proton  spectrum  which  is  quite  poor  in  all  areas.  The  combined 
chi-squared  for  this  calculation  increased  by  15%  over  the  best  fit 
value.  As  in  the  case  of  the  first  calculation,  the  0/A1  ratio  is  con¬ 
sistent  with  Al^O^  stoichiometry  but  the  depth  of  the  oxide  layer  is 

2 

about  30  yg/cm  greater.  This  increased  depth  with  no  hydrogen  is 
explained  wholly  by  the  hydrogen  contribution  to  the  stopping  cross 
section. 

Besides  giving  the  generally  good  detail  of  the  oxide  layer  and 
interface  as  described  above,  this  analysis  can  also  be  used  in  deter¬ 
mining  the  oxide  layer  density.  However,  such  surface  layer  density 
determinations  require  separate  measurement  of  the  layer  thicknesses 
by  techniques  such  as  optical  or  scanning  electron  microscopy,  or 
electrochemical  methods.  The  layer  density  is  then  given  by  p=d/t, 

2 

where  d  is  the  layer  areal  density  in  g/cm^  obtained  by  the  ion  beam 
methods  described  in  this  report  and  t  is  the  layer  thickness  in  cm. 

This  type  of  density  determination  is  currently  being  pursued  for  this 
anodized  aluminum  sample. 

The  concentration  profiles  shown  in  figures  2,  4  and  6  are  actually 
stepwise  functions  in  the  numerical  calculations.  The  step  size  varies 
over  the  calculation  depth,  being  small  where  any  concentration  gradient 
is  large  and  being  large  where  all  concentration  gradients  are  small. 

For  clarity  of  presentation,  these  stepwise  functions  are  represented 
by  smooth  curves  in  all  figures  showing  concentration  profiles. 

At  this  point,  the  question  of  uniqueness  of  the  fits  and  absolute 
accuracies  of  the  results  will  be  discussed.  As  was  pointed  out,  equally 
good  fits  to  the  alpha  spectra  of  figure  7  could  be  obtained  with  dif¬ 
ferent  concentration  profiles.  This  aspect  of  the  analysis  has  been 
noticed  with  many  other  spectra.  For  this  reason,  it  is  necessary  to 
fit  more  than  one  set  of  data  with  the  same  profile  set.  Thus,  requiring 
the  C,N,0  peaks  in  the  case  of  figures  3  and  5  and  a  fit  to  the  proton 
spectrum  in  figure  7  adds  considerably  to  the  reliability  of  the  re¬ 
sults.  To  definitively  test  the  uniqueness  of  a  solution,  a  known 
standard  sample  should  be  analyzed,  but  this  has  not  yet  been  done. 

The  absolute  accuracies  of  both  the  concentrations  and  layer 
depths  depend  totally  on  the  accuracy  of  the  cross  sections  and  stopping 
powers.  Also  as  the  layer  thicknesses  increase,  uncertaintites  in  the 
theoretical  expressions  for  energy  loss  straggling  become  more  important. 
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especially  in  determining  the  interface  profiles.  The  accuracies  of 
stopping  powers  are  presently  in  the  range  5-10%  so  the  layer  thickness 
determination  has  a  similar  uncertainty.  Typically,  elastic  cross 
sections  have  absolute  uncertainties  in  the  order  2-10%  while  reaction 
cross  section  uncertainties  range  from  10-20%.  The  statistics  in  the 
data  also  add  to  the  error  in  the  determination  of  concentrations.  In 
the  cases  of  the  C,N  and  0  reaction  peaks,  the  statistics  can  produce 
50%  uncertainties  in  the  concentration  profiles  in  some  regions.  How¬ 
ever,  since  those  concentrations  are  in  the  range  of  a  few  atomic  per¬ 
cent,  those  50%  uncertainties  only  amount  to  a  few  atomic  percent 
error  in  the  profiles.  This  result  is  generally  borne  out  since  the 
profiles  usually  need  to  be  adjusted  to  within  these  limits  for  the 
best  fits.  In  the  cases  where  elastic  distributions  are  fit,  the  cross 
section  uncertainties  dominate  so  that  the  errors  are  in  the  order  of 
2-10%  or  again  a  few  atomic  percent  in  the  concentration  profiles. 


IV.  SUMMARY 

A  general  and  versatile  program  which  was  developed  for  obtaining 
concentration  profiles  of  constituents  of  surface  layers  has  been 
described.  Its  usefulness  in  defining,  in  great  detail,  both  surface 
as  well  as  interface  structures  has  been  demonstrated.  Use  of  this 
analysis  program  on  spectra  generated  by  ion  beam  bombardments  will  aid 
in  answering  many  detailed  questions  abour  surface  layer  stoichiometries, 
buildup  of  oxide,  nitride  or  carbide  layers,  the  extent  of  carburiza¬ 
tion  of  steel  and  in  conjunction  with  optical  measurements,  the  density 
of  surface  layers. 
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APPENDIX  A 


-  PROGRAM  LISTING 
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C  PROFILE  —  PROGRAM  TO  ANALYZE  ROTH  RUTHERFORD  8ACKSCATTER I NG 

C  AND  NUCLEAR  REACTION  SPECTRA.  USES  AVAILABLE  VALUES  OP  CROSS 

C  SECTIONS  AND  ENERGY  LOSSES  AND  TAKES  ENFRGY  RESOLUTIONS  AND 

C  STRAGGLING  INTO  ACCOUNT.  THE  END  RESULTS  ARE  DEPTH  CONCENTRATION 

C  PROFILES  OF  CONSTITUENTS  OF  THE  SURFACE  OR  FILM  IN  ATOMIC  PERCENT. 

C  HISTORY  OF  THIS  CODE  - 

C  1974  -  ORIGINAL  WRITTEN  BY  J.E. YOUNGBLOOD  FOR  SEL-86? 

C  1975-76  -  MODIFIED  BY  A.NIILER. 

C  1977  -  MODIFIED  FOR  CDC7600  BY  J.GERRITS. 

C  1978  TO  NOW  -  CONTINUING  MODIFICATION  BY  A.NIILER. 

C  SUBROUTINE  INFORMATION. 

C  PROFILE  -  MAIN  ROUTINE  WHICH  DIRECTS  THE  FLOW  OF  TRAFFIC.  SFTS 

C  UP  GEOMETRIES  AND  HAS  MOST  OF  THE  INPUTS. 

C  SETTaB  -  INPUTS  AND  SETS  UP  CONCENTRATIONS.  CAN  OUTPUT  ENERGY 

C  LOSS  TABLES  FOR  CHECKING  PURPOSES. 

C  SET  I  ON  -  CALCULATES  IONIZATIONS  POTENTIALS 

C  WIDTH  -  CALCULATES  ENERGY  LOSS  STRAGGLING 

C  PICKE  -  PICKS  ENERGY  VALUES  TO  BE  USED  FROM  GAUSSIANS 

C  FINOE  -  NEEDED  BY  PICKE 

C  DIST  -  DISTRIBUTES  CALCULATED  YIELD  INTO  FINAL  ENERGY  DISTRIBUTION 

C  KIN  -  TWO  BODY  KINEMATICS  CALCULATION 

C  SIG  -  PEACTTON  OR  RUTHERFORD  CROSS  SECTIONS  ARE  DETERMINED 

C  FROM  TABLES  OR  CALCULATIONS 

C  SHAPE  -  DETERMINES  THE  ENERGY  AND  ENERGY  SPREAD  OF  PARTICLES 

C  TRAVELING  THROUGH  SURFACE  LAYERS 

C  FLOSS  -  CALCULATES  ENERGY  LOSEES  OF  HYDROGEN  AND  HELIUM  ISOTOPES 

C  MODIFICATION  NECDSSARY  TO  HANDLE  OTHER  TYPES  OF  PARTICLDS 

C  OUTPUT  -  EXPERIMENTAL  DATA  INPUT.  LEAST  SQUARED  FIT  TO  CALCULA- 

C  TED  DISTRIBUTION.  BOTH  ARRAYS  PREPARED  FOR  OUTPUT 

C  IGAUS  -  USED  BY  DIST 

C  PLOT  1  -  LINE  PRINTER  PLOT  OF  EXPERIMENTAL  ANO  CALCULATED  DISTS 

C  SWITCH  -  OUTPUTS  SWITCH  INFORMATION 

C  CONCTAB  -  OUTPUT  TABLE  OF  LAYER  CONCENTRATIONS 

CS  DEBUG 

C*  ARRAYS 

PROPPAM  PROFILE  (INPUT. OUTPUT .  X.APE5=  I ND  JT  ,TAPE6=0UTPUT» 
*DEBUG=OUTPUT) 

REAL  M1.MB.M3.M4 

C0MM0N/REACT/M1.M?.M3.M4.Q,ZIN(2) .ZTAR3(5) 

COMM ON/TABLES/ AC 0(5. 100.2) , CONCEN ( 5 ♦ 1 0 0 ) .CATARG(5.100).DFAC(100) 
COMMON  /FLAGS/FLAG.FLAGA.FLAGB.FLAGC. PLAGD.FLAGE 
COMMON/SHA/NSLA (100) .WT(IOO) 

COMMON/OU/X (410) ♦ Y ( 4 1 0 ) . Y Y ( 41 0 ) .XX(410) 

COMMON/S/  VALUE(2000) .ISSW(IO) 

C0MM0N/TIT/TITLE(20) .NP 
COMMON/ OEX/EM IN. EMAX 

DIMENSION  THICK (I  00) »F ION (100) .THKVAR(IOO) 

DIMENSION  SMI (5) .SM2(5) .SM3(5) .SM4(5) . SQ ( 5 ) . SE X ( 5 ) 

DIMENSION  THKIN(IOO) .THKOUT(IOO) 

INTEGER  ZIN.ZTARG 
DATA  FLAGCS/O.O/ 

C 

1  FORMAT ( 8F 1 0 . 0 ) 

2  FORMAT ( 10 15 ) 

3  FORMAT (20A4) 

II  FORMAT (IH1.25X.46HCALCULATI0N  OF  DEPTH  PROFILE  FROM  THE  SHAPE  OF, 
*3 1 H  THE  REACTION  PRODUCT  SPECTRUM  //// 
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*15H  REACTION  DATA  .10X.20A4/ 

*1  OX ♦ 4HM1  e,F10.6, 10X.4HM2  ■ , F 1 0 . 6 » 1  OX , 4HM3  =.F10,6.10X, 

*4HH4  *♦  FI  0  #  6/ 

*10X,19HREACTION  Q-VALUE  *  *  FI  0 . 6. 1 0 X  *  1 9HEXC IT  AT  I  ON  ENERGY  =,F10.5. 
*10X» 13HNET  Q-VALUE  *,F10.5/ 

*1  OX ,  1 3HBEAM  ENERGY  *»F10.5»22H  MEV  WITH  A  SPREAD  OF  ,F10.5,4HKEV  / 
*  10X, 14HTARGET  ANGLE  * » FI  0 . 3, 15X . 1 8HSC A TTER I NG  ANGLE  «,F10.3/ 

*1  OX , 53HTHE  REACTION  PRODUCT  DISTRIBUTION  WILL  BE  STORED  FROM, 

*F1 0 , 5  * 6HME  V  TO*F10.5*3HMEV// 

*15H  GEOMETRY  DATA// 

*1  OX ♦ 3BHTHE  BEAM-LIMITING  SLIT  HAS  A  RADIUS  0F.F10.5, 

IS.F10.5.20H  CM,  FROM  THE  TARGET/ 

DEFINING  SLIT  IS  REPRESENTED  BY,I5,7H  POINTS/ 

RADIUS  OF  THE  BEAM  IS,F10.5,18H  CM,  AT  THE  TARGET/ 
TARGET  IS  REPRESENTED  BY,I5.7H  POINTS/ 

DETECTOR  APERTURE  I S . F 1 0,=  . 1 1 H  CM.  AND  IS.F10.5, 

THE  TARGET) 


*1  OH  CM, AND 

*  1 5X , 35HTHE 

*  1  OX ♦ P5HTHE 

*  1 5X ♦ 28HTHE 
*1  OX » 24HTHE 
*?0H  CM,  FPOM 


1? 


FORMAT (/♦ 

*15X,40HTHE  DETECTOR  APERTURE  IS  REPRESENTED  RY=,I5,7H  POINTS/ 
*10X,?2HDETECTOR  RESOLUTION  IS,F10.5,4H  KFV// 

*  1 3H  TARGET  DATA  // 

COMPOSITION  OF  THE  TARGET  IS  DESCRIBED  BY, 15. 

OF  DIFFERENT  COMPOSITION/ 

FIRST  LAYER  IS  SUBDIVIDED  INTO  ,I5,12HEQUAL  SHEFTS/ 
DISTRIBUTION  OF  INCOMING  POINTS  AT  THE  CENTER  OF  FACH  , 
REPRESENTED  BY,I5.19H  OJTSOING  PARTICLES// 


*1  OX  ♦  45HTHF 
*32H  LAYERS 
*1  OX  #  35HTHF 
*  1  OX ♦ 58HTHE 
*?3HSHEET  IS 


*34H  THE  COMPOSITION  OF  THF  TARGET  IS  // 
*15X.10HIONIZATION»5(20H  ATOMIC) / 

*5X.SHLAYER*5Xf 11H  POTENTIAL  »5(20H  FLEM 


PERCENT) /) 


171  FORMAT (/.5X.26HIMP0SSIBLE  TARGET  LOC AT  I  ON , 2 X , 1 6H  ERROR  GEOM  FXIT) 
509  FORMAT ( 1 H 1 ) 

511  FORMAT (12X»47H***  PARTICLE  ENERGY  S°ECTRUM  PARAMETERS  *<»*.21X, 
123H**  DFPTH  PROFILE  **// 

*6X»3HED0»6X»5HSUM-S»6X»3HEDX»7X»3HKIN»5X»6HSUM-SP» 

*7X,1wE»6X,6HSIGEDX,2X»4HRH0X,7X,5H0EPTH,12X,3HYTE/ 

*5X  » 5H  ( ME V  )  *  5X  » 5H  ( MEV )  »5X.5H(MEV)  » 1 5X  *5H  ( ME V  )  « 5X  ♦  5H  ( ME V  )  *  3 X  ,  7H  <  mr/s 
*R ) *  7X  *  BH ( MG/CM2) . 3X . 9H ( MICRONS) , 7X, 0R ( AT/CM? ) ) 

51?  FORMAT  (f.F10,4.F10,3*4X.F10,4»F13,5.El  5.5) 

513  FORMAT ( 13HREAR  OF  L AYER . 1 4 ♦ F 1 1 . 4 . FI 3 .4 ♦ F 1 0 . 4 . F44 . 5 . F 1 3 . 2, F 1 7 . 5/ 
*50X*37HOUTPOING  PARTICLE  CANNOT  LEAVE  TARGET/) 

FORMAT (9X.5HL AYER* I4,F78.5,F13.2,F17.5) 

FORMAT ( 1H1 ) 

FORMAT (1HO»38HREACTION  IS  DOUBLE- VALUED  OR  FORBIDDEN) 

FORMAT ( 44H  INCOMING  OR  OUTGOING  PARTICLE  GRAZES  TARGET) 

FORMAT ( 29H  DETECTOR  IS  SHADOWED  BY  SLIT) 

FORMAT ( 1H0.50HFNERGY  LOSS  STRAGGLING  CALCULATION  IS  NOT 


514 

523 

526 

527 

528 

529 
613 


1000 

1111 

4771 

C 


RELIABLE) 

FORMAT (39HT0TAL  CONCENTRATION  OF  TARGET  COMPONENT , F8. 2* 1 H= . El b . 5 • 
♦7HAT/CN2. ) 

FORMAT (315. 5F 12,5. IX, 2 (1PE12.5) ) 

F  ORMAT ( 8E 1 2  »  5 ) 

FORMAT ( BH  THKVAR ( . I5.4H)  c  ,F10,5.5H  STOP) 


9991  READ(5»3) TITLE 
DO  141  1=1,100 
141  THKVAR(I)bO. 

FL AGCS* 1.0 

DTOR=3. 141 592654/130. 
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ABO  =  80.0*DTOR 

INPUT  DATA  FOR  CALCULATION  OF  CONSTANT 

BEAM  CHARGE  IS  GIVEN  BY  CHAR.  DETECTOR  SOLID  ANGLE 

GIVEN  BY  OMEGA. 

NCONSsO 
1  =  1 

CALL  SWITCH 

READ ( 5* 1 )  CHAR*  OMEGA 

10  READ (5* 1) M1,M2,M3»M4»QRE ACT, EXCIT 
IF(Ml.NE.0.)GO  TO  414 

Ml* SM 1(1-1) SM2*SM2 ( 1-1 ) $M3«SM3 ( 1-1 ) $M4=SM4 ( 1-1 ) 

QREACT=SQ ( 1-1 ) SEXCIT=SEX( 1-1) S60T0411 
14  SM1(I)*M1 

SM?(I)=M2 
SM3 ( I ) =M3 
S*4 ( I ) =M4 
SQ  ( I ) =QREACT 
SEX ( I ) =EXCIT 
Q=SO ( I ) -SEX ( I ) 

1*1*1 

NC0NS=NC0NS  +1 
GO  TO  410 

11  CONTINUE 

M1,m?,M3.M4.QREACT  DEFINE  THE  NUCLEAR  REACTION 
E  IS  THE  ENERGY  OF  THE  INCIDENT  BEAM  ( ME V ) . 

THETAT  IS  THE  ANGLE  BETWEEN  BEAM  AND  TARGET  NORMAL. 

THETAD  IS  THE  SCATTERING  ANGLE  FOR  THE  DETECTOR  CENTER. 

EXCIT  IS  THE  EXCITATION  ENFRGY  (MEV). 

EMAX  AND  EMIN  DEFINE  THE  LIMITING  ENERGIES  OF  THE  0UT°UT  ARRAY. 
DTHETA  ALLOWS  ANGULAR  SPRFAO  OF  OUTGOING  PARTICLE  TO  BE  FOLDED 
INTO  A  CALCULATION 

READ (5* 1 ) E. THETAT* THETAD. EMAX, EMIN. DTHET 
RS  AND  DS  ARE  USED  TO  DEFINE  DIVERGENCE  OF  REAM.  IN  THIS  CASE, 
THEY  ARE  TAKEN  AS  A  fully  illuminatfd  quaorupole  aperture. 

RT  IS  THE  TARGET  RADIUS  AT  NORMAL  INCIDENCE.  IN  THIS  CASE. 

RT  IS  THE  RADIUS  OF  THE  ANTI-SCATTFR ING  SLIT. 

RD  AND  DD  DEFINE  THE  DETECTOR  LOCATION. 

READ  (5,1)  RS , DS , RT , RD , DD 

NGS  IS  THE  NUMBER  OF  POINTS  USED  TO  DEFINE  THE  SLIT. 

NGT  IS  THE  NUMBER  OF  POINTS  USED  TO  DEFINE  THE  TARGET. 

NGD  IS  THE  NUMBER  OF  POINTS  USED  TO  DEFINE  THE  DETECTOR . 

READ(5,?)NGS,NGT.N6D 

DETRES  IS  THE  DETECTOR  RESOLUTION  IN  <E V 
ERES  IS  THE  BEAM  ENERGY  RESOLUTION  IN  KEV 
PEA0(5,1)DETRES,ERES 

NLAYER  IS  THE  NUMBER  OF  LAYERS  IN  THE  TARGET.  LIMITED  TO  100. 
NSLAR  IS  THE  NUMBER  OF  SECTIONS  OF  THE  FIRST  TARGET  LAYER.  ALL 
SUBSEQUENT  LAYERS  WILL  HAVE  SLABS  0?  THE  SAME  THICKNESS. 

THIS  IS  DONE  IN  4730  LOOP. 

NDIV  IS  THE  NUMBER  OF  POINTS  USED  TO  REPRESENT  THE  GAUSSIAN  SHAPF 

READ (5, 2) NLAYER, NSLAB.NDIV 

IF(NLAYER.GT.IOO)  nlayer=ioo 

IF  (NSLAR. GE. 100)  NSLA8=100 

ANGO=THETAD*DTOR 

ANGTxTHETAT*DTOR 

I F ( ISSW (4) .NE.l )  GO  TO  5 

WRITE (6, 1 1 ) TITLE, Ml , M2, M3.M4.QREACT, EXCIT, Q,E, ERES. THETAT, THETAD, 
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♦EMIN.EMAX.RS.DS.NGS.RT.NGTiRD.DD 
WRITE (6. 12) NGD.DETRES. NLAYER. NSLAB.NDIV 
5  CONTINUE 

INPUT  THE  CONCENTRATION  DATA  AND  GET  ENERGY  LOSS  INFORMATION 
OUTPUT  IF  ISSW<1)*1 
CALL  SETTAB(NLAYER.THICK.FION) 

PROGRAM  RETURNS  HERE  AT  COMPLETION  OF  OUTFUT. 

INITIALIZE  THE  OUTPUT  ARRAY  'VALUE'. 

CALL  OUT  1 

THICK  IS  THE  ACTUAL  THICKNESS  OF  EACH  LAYER  IN  MG/CM**2 

USE  A  NEGATIVE  VALUE  OF  THICKilJ  TO  INITIATE  OUTPUT  OF  DATA. 

ONLY  ONE  CARD  NEED  BE  USED  TO  FLAG  OUTPUT  EVEN  IF  NLAYER 
EXCEEDS  8. 

READ ( 5.  *)  ( THICK ( I ) *I*1.NLAYERI 

*  *  *  *  »•••■*»»*«» 

DO  GO  1*1 .NLAYER 
30  THKVAR(I) *0.0 

NO  THICKNESS  VARIATION  IS  ALLOWED  AS  FZ  12DEC78.  THESE  FEW 
CARDS  NEED  TO  BE  UPDATED  SHOULD  TOE  NEED  ARISE  TO  START  AGAIN 
WITH  SURFACE  ROUGHNESS  CALCULATIONS) 

*  *  *  *  ##**■#*###* 

DO  4730  ISLA*1 .NLAYER 
II*  ISLA  *1 

IF  (  THKVAR  (II)  .NE.0.)WRITE(6*4T71)  II»TR<VAF:(  II) 

4730  NSLA ( ISLA) *THICK( ISLA) /THICK ( 1) *NSLAB 
IF ( THKV AR ( 1 ) . NE . 0 . ) NSLA ( 1 ) * 1 
NSL AA*NSLAB 

IF ( ISSW (6) ,NE.0)WRITE(6.523) 

BEGIN  COMPUTATION 

FLAGxO.O 

FLAGA*0.0 

FLAGB*0.0 

FLAGC*0.0 

FLAG  NOT  0  INDICATES  THE  RF  ACT  I ON  CANNCT  OCCUR  OR  IS  OOUBLE-V ALUFD 
FLAGA  NOT  0  INDICATES  THAT  INCOMING  OR  OUTGOING  PARTICLE  GRAZES  TARGET 
FLAGS  NOT  0  INDICATES  SHAOOWING  OF  DETECTOR  RY  SLIT 

FLAGC  NOT  0  INDICATES  THE  GAUSSIAN  WIDTH  IS  NOT  RELIABLY  CALCULATED 
NIX»2.*DTHET*1 
DTHET«DTHET*DTOR 
ANGT  *  ANGT  -DT MET 
DO  600  MULT  *  l.NIX 

xmult*mult-i 

ANGT  *  ANGT  ♦  XMULT  *  DTOR 
DO  500  IGS*1 . NGS 

C  SELECT  A  POINT  IN  THE  SLIT. 

XS*DS 

YS*0. 

IF ( NGS.LT •?)  GO  TO  51 

YS«RS*(1.-2.*FLOAT(IGS-1)/FLOAT(NGS-1) > 

51  CONTINUE 

DO  500  IGT*1 »NGT 

C  SELECT  A  POINT  IN  THE  TARGET. 
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XT*YT*0  • 

IF ( NGT »LT.?)  GO  TO  52 

XT*RT*SIN(ANGT)*(2.*FL0AT(IGT-1) /FLOAT (NGT-l)-l.) 

YT«RT*COS( ANGT)* (1 .-2. *FLOAT ( I GT-1 ) /FLOAT < NGT- 1 ) ) 

C  THE  NEXT  2  STATEMENTS  ASSUME  RT  IS  * DEr INED •  BY  BEAM  SIZE. 

XT* XT /COS (ANGT) 

YT*YT/COS(ANGT) 

52  CONTINUE 

DO  500  IGD*1 * NGD 

C  SELECT  A  POINT  IN  THE  DETFCTOR  APFRTURE. 

XD*-DD*COS< ANGD) 

YD*DD*SIN(ANGD) 

I F ( NGD # LT * 2 )  GO  TO  53 

XD=-OD*COS< ANGD) ♦ RD*S IN ( ANGD) *< 2. "FLOAT (IGD-1 ) /FLOa T ( NGD- 1 ) -1 . ) 
YD*DD*SIN(ANGD) +RD*COS ( ANGD )*( 2 . "FLOAT ( IGD- 1 ) /FLOAT ( NGD- 1 ) -1 .) 

53  CONTINUE 

IF ( XS.LE. XT)  GO  TO  170 
YDAsYT* (RS-YT) * (XO-XT) / (DS-XT) 

IF (YD.LE.YDA)  GO  TO  180 
GEOA.(XT-XS) **2+ ( YT-YS)"*2 
GEORs(XT-XD) **2* (YT-YD) *"2 
GEOC* (XO-XS) *"2  + ( YO-YS) **2 
U=GEOA+GEOP-GEOC 
U=U"0.5/SQRT(GEOA*GEOB) 

ANGLAB*AC0S(AMIN1 <U,1  .)  ) 

ANGLAB«3.14159-ANGLAB 
ANGIN*ANGT-ATAN( ( YS-YT) / (XS-XT) ) 

ANG0UT*3.141592654-ANGLAB-ANGIN 

deglab*anglab/dtor 
degin*angin/dtor 

DEGOUT *AN GOUT /OTOR 

GEOMFTRIC  PARAMETERS  HAVE  BEEN  DEFINED  FOR  IGS 1 1 GT  »  AND  IGD. 

TEST  FOR  VALID  GEOMETRY 

IF  ( ANGIN.LT.A80. AND. ANGOUT.LT. A80)  GO  TO  100 

PARTICLE  FNTERS  OR  LEAVES  WITHIN  10  DEGREES  OF  SURFACE  PLANE. 
FLAG  GRAZING  INCIDENCE  AND  SKIP. 

FLAGA*1 .0 
GO  TO  500 
170  WR I TF ( 6 » 1 7 1 ) 

BAD  GEOMETRY  FOUND. 

STOP 

180  CONTINUE 

DETECTOR  SHADOWEO  BY  SLIT.  FLAG  AND  SKIP. 

FL AGR* 1,0 
GO  TO  500 

GEOMETRY  IS  VALID. 

100  CONTINUE 

U=1.0/COS(ANGIN) 

V=1.0/COS( ANGOUT ) 

DO  110  1*1 tNLAYER 
THKIN(I)*THICK(I)"U 
110  THKOUT (I)*THICK(I)"V 

IF (THKVAR(l) . EO.O. )  GO  TO  7878 
THKIN(1)*0.5"THKIN(1) 
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THKOUT (1)»0.5*THKOUT (1) 

7878  CONTINUE 
WIDT.O. 

DO  498  KM*1 * NCONS 
Ml«SMl (KM) 

M?*SM2 ( KM) 

M3«SM3(KM) 

M4*SM4(KM) 

IF(M2.EQ.0> GO  TO  498 
Q»SQ(KM)-SEX(KM> 

DO  499  LLAYER»1 *NLAYER 

layer»llayer 

IF (CONCEN(KM» LAYER) .LE.O.) O0T0499 

EA  is  the  mean  energy  and  wid  is  the  standard  deviation  OF  THE 

INCIDENT  BEAM  ENERGY  DISTRIBUTION. 

INITIALIZE  EA  AND  WIDSQ  FOR  THIS  LAYER. 

CALL  SHAPE (E* DE* WID* Ml ♦ l.F ION* THKIN.THK VAR t LAYER) 

IF(E.LE.DE) GOT  0499 
EAsE-DE 
EAO*EA 
W I DS A»0 , 

EAO  IS  THE  ENERGY  OF  THE  INCIDENT  PARTICLE  AS  IT  ENTERS  THE  L A YFR 

W  IDSQ*W  ID**2  ♦  WIDT 

THKBsO.O 

UAr  1  , 

NSLARaN SLA (LAYER) 

THKN.THKIN (LAYER) /NSLAB 

EACH  LAYER  OF  INTEREST  IS  DIVIDED  INTO  NSLAB  SLABS. 

IF (LAYER. EO.l.AND.THKVAR(LAYER) . NE. 0 , ) SO  TO  4177 
GO  TO  4178 

SKIP  TO  ST  4178  IF  NO  SURFACE  ROUGHNESS  IS  USED. 

177  LAYER-2 

************* 

THIS  REGION  TO  BE  USED  IF  SURFACE  ROUGHNESS  CALCULATION  DONE  AGAIN 
«*»»*»»«»»*»« 

4178  CONTINUE 

DO  300  NS«1. NSLAB 
THK  A»THKN*U A 
THKRrTHKR*THKA 

CALL  EL0SS(EA*LAYER*1.ThKA.DE) 

IF ( E A, LE , DE ) G0T0499 

CALL  WIDTH(EA. DE. Ml ,FI0N (LAYER) .WIDA. LAYER) 

WIDSA»WIDSA*WIDA**2 

FA«EA-DE 

WIDsWIDSA  ♦  WIDSQ  ♦ ( ERES* , 00 1 > **2 

THE  INCIDENT  BEAM  NOW  HAS  A  GAUSSIAN  01 STR I  HUT  I  ON ( E A . W I D ) . 
WIDsSQRT(WID) 

THE  INCIDENT  PARTICLE  HAS  REACHED  THE  REACTION  SITE. 

REPRESENT  THE  DISTRIBUTION  BY  SELECTING  NDIV  ENERGIES  SUCH  THAT 
THE  INTEGRAL  OF  THE  DISTRIBUTION  FROM  ONE  ENERGY  TO  THE  NEXT 
YIELDS  A  CONSTANT  AREA. 

DO  200  ND*1 »  NDI V 

UB= (FLOAT (ND) -0.5) /FLOAT (NDIV) 

UB*1 . 0-UB 

C  SELECT  NDIV  REPRESENTATIVE  ENERGIES  SJCH  THAT  THE  INTEGRAL 
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C  OF  THE  DISTRIBUTION  from  -INFINITY  TO  E  EQUALS  1  -  ( N- 1 /? ) /ND I V . 

CALL  PICKE  (UB.EAfWID.EB) 

C  THE  REACTION  NOW  OCCURS.  THE  OUTGOINS  PARTICLE  HAS  AN  ENERGY  EC. 

CALL  KIN(EB.ANGLAB.EC.XSEC.KM) 

IF ( XSEC.LE *0.0)  00  TO  300 

THE  DISTRIBUTION  OF  THE  OUTGOING  PARTICLE  IS  NOW  DETERMINED  AS  IT 
TRAVERSES  THE  sample. 

THKC*THKA*THKOUT (LAYER) /THK IN (LAYER) 

WIDSQB«0. 

DO  402  KK»1,NS 

CALL  ELOSS(EC.LAYER.2.THKC.DE) 

IF(EC.LE.DE)  GO  TO  300 

CALL  WIDTH(EC.DE.M3,FI0N(LAYER) .WIDB. lAYER) 

WIDSGB«WIDSQB*WIDB**2 
EC=PC-DE 
40?  CONTINUE 

THE  REACTION  PRODUCT  HAS  LEFT  THE  LAYER  IN  WHICH  THE  REACTION 
OCCURRED. 

ED«EC 

CALL  SHAPE(ED.DE. WIDB. M3,?, FI  ON. THKOUT.THKVAR, LAYER) 

IF (ED.LE.DE)  60  TO  300 
WIDSQB»WIDSQB*WIDB**2 
IF (NDIV.EQ.I) WIDSQB=WIDSQR+WID**2 
INCREASE  WIDB  TO  INCLUDE  THE  EFFECT  0^  DETECTOR  RESOLUTION. 
WlDBrSQRT(WIDSQB+ ( 0 . 0 0 1 *DETRES ) **2) 

EE=ED-DE 
GO  TO  3178 

««««««««««««« 
AGAIN.  STUFF  RELATED  TO  SURFACE  ROUGHNESS  REMOVED 
####**#*#  »  »  « 

3178  CONS«CHAR*OMEGA*3.906» (10. 0**6) 

WEIGHT *XSEC  *  CONCEN ( KM, LAYER) *TH« N*CONS 
THKN  IS  THE  SLANT  RANGE  THROUGH  THE  SLAB  FOR  THE  INCIDENT  BEAM. 
THE  DISTRIBUTION  OF  THE  OUTGOING  PARTICLE  AT  THE  DETECTOR  IS 
GAUSSIAN  (EE.wIDB). 

THE  RELATIVE  WEIGHT  OF  THIS  DISTRIBUTION  IS  'WEIGHT*. 

THE  DISTRIBUTIONS  ARE  ACCUMULATED  IN  THE  OUTPUT  ARRAY ' VALUE • . 

CALL  DIST(EE.WIDB.WEIGHT.EMAX.EMIN) 

IF(ISSW(6) .NE.l. OR.EE.LT. EMINJGO  TO  200 
WRITE  (6. 1000) LAYER. NS, ND. EA.WID.EE, WIDB. 

S  CONCEN(KM, LAYER) .WEIGHT, XSEC 
00  CONTINUE 

00  CONTINUE 

99  CONTINUE 
98  CONTINUE 
00  CONTINUE 
600  CONTINUE 

I F ( ISSW (10) • EQ. 1 ) CALL  CONCTAB ( NLAYER. THICK) 

NSLABaNSLA (1) 

NP*0 

CALL  OUTPUT ( YMX ) 

CALL  0UT2 

THE  PUBLIC  RELATIONS  BLURB  HAS  BEEN  REMOVED  DUE  TO  RELATIVELY 
NONSENSICAL  INFORMATION  CONTAINED  THEREIN. 

OUTPUT  FLAGS 
WRI TE ( 6 . 523 ) 

IF (FLAG.GT.O. 1)  WRITE(6,5?6) 
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IF(FLAGA.GT.O.l)  WRITE(6,527) 
IFIFLAGB.GT.O.J)  HRITE (6*528) 
IF(FLAGC.GT.O.l)  WRITE<6.529> 
RE AD ( 5*  2) I 

I F  d . NE .0 ) GO  TO  9991 

STOP 

END 


SUBROUTINE  SETION (Z.FION.N) 

THIS  SUBROUTINE  CALCULATES  THE  EFFECTIVE  IONIZATION  POTENTIAL 

INTEGER  Z 
REAL  ION 

COMMDN/SET/ATARG(S.100) . ATMASS( 100) 

DIMENSION  IONUOO)*Z<1) 

DATA ( ION ( I )» 1*1 *94) /1 8*  5. 42* 0 *37. 8*66. 8 » 73. 0 »  77 .5 » 84*  7*  88 • 8 *  1 21  •  0  » 
*131. 0*1 43. 0.152. 4. 162 .5. 172. 2.181*5.1 92. 0.200. 6. 208. 8*218. 5. 228.0. 
*237. 3 *246. 4 *255. 3 *265. 2 *275. 0*284. 7 *294. 3 *304. 1.313.8*324.0. 333.9. 
*343.7*353.1 .362. 1.371. 0.380 .9. 390. 7 .400. 5 .4 10. 3. 420. 0.429. 7. 439. 3. 
*448.9.458.5*  468 .0.477 .9.487.8.497.7 .507 .6.517 .5.527 .3.537 .2.547 .0 . 

*556,7.566,5.576.2.586.0.595.7*605.3.61 5.0.624.6.634.3.643.9.653,4, 
*663.0.672,9.682.7.692.6,702,4,712.3.722. 1 .731.9.741 .7.751 .5.761 .3. 
*771. 0*780.8*790.5* 800. 3*810. 0*819.7*829.4, 839, 1,848.8, 858. 5, 868. 2, 
*877.8,887.5,897.1 .906.8,916 .4. 926. C .935.6,945./ 

U=V*W*0 . 0 
DO  10  1*1.5 

IF(ATARG(I,N) .LT.O. 0001)60  TO  10 
J*Z(I) 

THE  DENSITY  IS  PROPORTIONAL  TO  Ad).  T HF  RELATIVE  WEIGHT  OF  A 
CONSTITUENT. 

X=FLOAT(J)*ATARG(I.N) /ATMASS ( J) 

U  =  U  ♦  x 

W  =  W*X*ALOGdON(J)  ) 

10  CONTINUE 
WsW/U 

F ION*EXP ( W) 

RETURN 

END 
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SUBROUTINE  SETT AB ( NL AYER , THICK » F ION) 

THIS  SUBROUTINE  INPUTS  AND  SETS  UP  THE  CONCENTRATION  TARLES  AND 
CALCULATES  AND  OUTPUTS  ENERGY  LOSS  TABLES  FOR  CHECKING  PURPOSES 
ISSW(l)  HAS  TO  BE  1  TO  GET  THIS  ELOSS  OUTPUT.  SEE  THF 
DESCRIPTION  IN  FRONT  OF  ELOSS  SUBROUTINE  FOR  HOW  ENERGY 
LOSSES  ARE  CALCULATED. 

COMMON/TABLES/ACO(5,lOO,2) » CONCEN ( 5 » 1 00 ) . CA TAR3 < 5 , 1 00 ) ,DFAC(100) 
COMMON/REACT/MI , M2 , M3 » M4 , Q, Z I N ( 2) »ZTAR3(5) 

COMMON/S/ VALUE (2000) • ISSW (10) 

COMMON/TSYMB/SYMBUOO)  .ICON 
C0MM0N/SET/ATARG(5. 100)  .ATMASSUOO) 

COMMON/ELO/IZA 

DIMENSION  THICK  ( 1  00  )  .FI0N(  100)  »EREF(37)  .DUMTUOO)  .  SI  (5.37.2) 
DIMENSION  S2(15) .ARA(IOOO) 

INTEGER  SYMB ♦ Z IN »  ZTARG 
PEAL  Ml, M2, M3. M4 

DATA  SYMB/2H  H. 2HHE . 2HL I . ?HBE . 2H  B.2H  C.2H  N.2H  0. 

12H  F.2HNE.2HNA.2HMG.2HAL.2HSI »2H  P.2H  S.2HCL.2H  A,2H  K . 2HCA , ?HSC » 
•2HTI.2H  V.2HCR.2HMN.2HFE.2HC0.2HNI.2HC0.2HZN.2HGA.2HGE.2HAS.2HSP, 
*2HBR.2HKR,?HRB.2HSR.2H  Y. 2HZR . ?HNB . 2HMD . 2HTC . 2HRU . 2HRH, 2HPD . ?H A G. 
*2HCD.2HIN.2HSN.2HSB.2HTE.2H  I , 2HXE , 2HCS , 2HB A . 2HL A , 2HCE . 2HPR , 2HN0 . 
#?HPM,2HSM.2HEU.2HGD.2HTB.?HDY.2HH0.2HER»2HTm,2hY8.2HLU.2HHF,?HTA, 
*2H  W.2HRE.2H0S.2HIR.2HPT.2HAU.2HHG.2HTL.?HPB.2HBI.2HP0.2HAT,?hRN. 
*2HFR.2HRA.?HAC.2HTH.2HPA.2H  U , 2HNP, 2H? J , 2HAH, 2HCM , 2H8K , 2HCF  , 

*2H  . 2H  / 


DATA 

DF AC/597. 7. 150,5.86,83.66,85 

,55. 

73*50. 

16,43.02,37.66 

.3 

1.71, 

*29.85 

.26. ?1 

.24,78.22.33 

•2lt45tl9* 

45,1 

9.79,1 

6.99,15.08,15. 

41 

,15.03, 

*13.4, 

12.58. 

11.83,11.59. 

10.97, 10.79, 10 

.22,10 

.26,9.482,9.21 

7, 

8. 642, 

*8.3, « 

.042.7 

.63,7.54,7.1 

9,7.049,6. 

876, 

6.777, 

6.605,6.485,6. 

28 

,6.066. 

*5.951 

♦5.855 

,5.663,5.585 

,5.36,5.247,5. 

076,4. 

949,4.722,4.746, 

4.569, 

*4.533 

.4.387 

,4,337.4,3.4 

.276.4.177 

,4.0 

99,4.0 

07,3.965,3.831 

,3 

.791 , 

*3 ,70R 

.3,653 

,3.602,3.567 

,3.482,3.443,3 

.376,3 

.33,3.272,3.236, 

3.168, 

*3.135 

.3.088 

,3.059.3.004 

,2.948,2.908,2 

.983,2 

.869,2.869,2.7 

14 

,2.702, 

*2.666 

.2.654 

,2.597,2.608 

,2.531, .0, 

.  0 ,  . 

0  ,  .  0 ,  . 

o*.o«.o«.o/ 

DATA  ( 

ARA(I) 

,1*1.190)/ 

S 

1. 

44, .45, .2426 

,12.0,115. 

9,1. 

397, .45, .4845, 5. R73, 

52 

.25, 

SI. 6,. 

45. ,7256.3.013.45, 

78,2.59, .45, .9 

66. .15 

38,34.75. 

$2,815 

..45.1 

.206,1.060,28.55,2.989 

.  .45 

,1.445 

..9572.28.19, 

$3.35. 

.45.1. 

683,1.90,25. 

13, 3.0, .45 

,1.9 

2,2.0, 

22.3, 

$2,352 

..45.2 

.157,2.634,1 

8.16,2.199 

,.45 

,2.393 

,2.699,15.68, 

$2,869 

..45.2 

.628,1.854,1 

4.72,4.293 

,.45 

,2.862 

.1.009,13.97, 

$4,739 

.  .45.2 

.766. .1645. 20. 23,4.7,. 

45,3 

.329, . 

55,13.21, 

$3,647 

. .45.3 

.561,1.56,12 

.67,3.691, 

.45, 

3.792, 

1.219,12.11, 

$5,714 

. . 45 . 4 

.023. .8786,1 

1.78,6.5,. 

45,4 

.253.  . 

53,11.23, 

$5,833 

. .45.4 

.482. .5457,1 

1.29,6.252 

.  .45 

,4.71, 

.5533,11.12, 

$5,884 

. . 45 . 4 

.938, .5609,9 

. 995 , 5.496 

.  .45 

,5.165 

,  .5695,9.474, 

$5,055 

. .45.5 

.391  ,.9523,9 

.117,4.489 

.  .45 

,5.616 

,1.336,8.413. 

$3,907 

. .45.5 

.725,1.461,8 

.829,3.963 

,.45 

,6.065 

,1.243,7.782. 

$3,535 

♦  .45.6 

.288,1.372,7 

.361 ,4.004 

,.45 

,6.205 

,.5551,8.763, 

$4,175 

. .45.4 

.673, .3878,2 

1.88,4.75, 

.45. 

6.953, 

.2952,6.809, 

$5,697 

. .45.7 

.173*. 2026. 6 

.725,6.3,. 

45,6 

.496,. 

11,9.689, 

$6.01? 

. .45.7 

.611, .2925, 6 

.447,6.656 

..45 

,7.395 

,.11 75 ,7 .684, 

$6,335 

..45*8 

.046. .3652,6 

. 244 ,7.25, 

.45, 

9.262, 

.22,6.087, 

$6,429 

. .45.8 

.478, .2929.6 

.087,7.159 

,.45 

,9.693 

.  .3303,6.003/ 

D  AT  A  ( 

ARA(I) 

.1*191,380)/ 

$7,234 

. .45*8 

.907. .3678,5 

.889,7.603 

,.45 

,9.12, 

.4052,5.765, 

$7,791 

..45*9 

.333. .4427,5 

.587,7.248 

37 

,.45 

,9.545 

.  .4902,5.367, 

57.671 ,  . 45 ,9.756 ♦.5176 » 5.315 » 6. 887 » .45 *9.966..  5551 ,5. 151 . 
*6.677, .45* JO. 18*. 5925* 4. 91*. 5. 9,. 45 ,10. 38 ,.63, 4. 758, 

56.354. . 45,7.160* .3376. 13. *4. 6. 554.. 45. 10.0. ,3555.4.626. 

S7. 024 ..45 .11, 01 ..3709.4,54.7, 227.. 45. 11. 21.. 3864. 4. 4  74. 
SR.4R, .45,8, 608, .348,9.074,7.871 ,.45.1 1.62,. 39 24 .4, 4 02. 

S8. 71 6*. 45, 11. 83, .3948,4.376,9.289, .45 ,12. 03*. 3973. 4. 384, 

58. 21 8. . 45,12.23, . 3997 » 4 .447 » 8 . 91 1 . . 45 , 1 2 . 43 » . 4021 , 4 . 51 1 . 
59.071, .45, 12. 63*. 4045, 4. 54, 8. 444, .45,12.83, .4069,4.42, 

58. 2 19..  45,J 3. 03,. 4093, 4. 298, 8.,. 45, 1 3 . 23, . 4 1 1 8, 4. 1 82, 

57.786. . 45. 13.43. .4142.4.058.7.58. .45. 13. 62.. 4 166. 3. 976, 

57.30. . 45. 13. 82.. 4 19. 3. 877. 7.591 . .45. 14. 02.. 42 14. 3. 863, 

56. 996..  45. 14. 21.. 4239. 3. 725. 6. 21.. 45. 14. 4 ..4263. 3. 632, 

55. 874..  45. 14. 6.. 4287. 3. 498. 5. 70 6. .45. 14.79. .433.3.467, 

55.542. . 45. 14.98. .4335. 3. 34 2. 5.386. .45. 15.17. .4359.3.292, 

55.505. . 45. 15.36. .4384.3.243.5.657. .45.15.55. .4408.3.195, 

55. 329..  45. 15. 74.. 4432. 3. 186. 5. 144. .45.15.93. .4424.3.144, 

55. 851..  45 .16. 12 ..44 16. 3. 122. 5. 704.. 45. 16. 3.. 440 9. 3. 002, 

DATA {ARA ( I) ,1*381,590)/ 

55. 563..  45. 16. 49.. 440 1.2. 965 .5. 034.. 45. 16. 67.. 4393. 2. 871, 

55.46. . 45. 18.32. .4385 . 2 .542 . 4. 843. . 45. 1 7 . 04 , . 4878 ♦ 2. 882 , 

55.311 .. 45. 17.22,. 537, 2.91 3, 5. 982, .45, 17.4,. 5R63, 2.871 , 

56. 7..  45 .17. 8.. 677. 2. 66. 6.928 ..45. 17. 7 7.. 5863. 2. 01 2, 

56.979. . 45. 17.95. .5863.2.776.6.954. .45. 10. 12. .5863.2.748, 

57.82. . 45. 18. 3.. 5863. 2. 737. 8.448.. 45. 18. 48. .5863.2.727, 

58. 6 09 .. 45. 18. 66.. 5863. 2. 697. 8. 679.. 45. 10. 83,. 5863, 2. 641  , 

58.336. . 45. 19.01. .5863.2.603.0.204. .45.19.18. .5863.2.573, 
50,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0*0,0,0,0,0,0,0,0,0,0,0,0,0,0, 
50,0,0,0,0,0,0,0,0,0, 

5.9661 .. 4126.6.92.8.831.2.582.2.027. .293 1.26. 34. 6. 66 ..3409, 

5 1.42..  49. 12. 25. 32. 9. 161.. 2. 206.. 51. 15. 32.. 25. 8. 995, 

53. 691..  4 120. 18. 48. 50. 72. 9.. 4. 232.. 3877. 22. 99. 35. .7.993, 

52.51 . . 4752.38.26. 13.02. 1 .892. 1.766. .5251 .37.1 1 .15.24.2.004, 

51. 533. . 53 1.40. 44. 18. 41. 2. 718.1. 183. .55. 39. 83 ,  17. 49, 4. 001, 

59.894. . 3001.23.65. .384. 92. 93. 4. 3,. 47, 3  A. 3, 3. 3, 12. 74, 

52.5. . 625.45.7. .1.4.359.2.1. .65.49.34.1 .788.4.133, 

51. 729. . 6562. 53.41. 2. 405.3. 045. 1.402.. 6791.58. 98. 3. 520.3. 211, 

51. 117..  7044 .69. 69. 3.705. 2. 156 ..91 72.. 724. 79. 44. 3. 648. 1.646, 
DAT  A ( ARA ( I ) ,1*591,780)/ 

58. 554.,  38 17. 83. 61. 11. 84, 1.875, 6. 297,. 4622. 65. 39, 10. 14,5.036, 

55. 307..  491 0.61. 74. 12, A, 6. 665, 4. 71 . . 50 97 ♦ 65. 2b , 8 . 806 . 5 . 940 , 

56.151 .. 4524.83. .18.31 .2. 71 .6. 57.. 4322. 84. 76. 15. 53.2. 779, 

55. 738..  4492. 84. 61 .14. 18.3.101.5.013. .4707.85.58.16.55.3.211 , 

54.32. . 4947 .76. 14. 10. 85. 5. 441. 4. 652,. 457 1,80. 73. 22. ,4.952. 

5  3.1 14 ,.5236 ,76.67, 7. 62, 6. 385, 3. 11 4, .5236,76.67,7.62,7.50  2, 

53.1 14..  5236. 76. 67. 7. 62. 8. 51 4. 5. 74 6. .4662.79.24. 1. 185. 7.993, 

52.792. . 6346.106.1 . .2986.2.331.4.667. .5095.124.3.2.102.1 .667, 

52. 44. . 634 6. 105. ..83. 2, 85 1,1. 491*. 7118, 120. 6, 1.101, 1.877, 

5 11. 72..  3826. 102. 0.9. 231. 4. 371. 7. 126.. 4 904. 11 9. 3. 5. 784. 2. 454, 

51 1.61..  3955. 146. 7. 7. 031 .1.4, 23, 10. 99.. 41. 163. 9, 7. 1,1. 0  52. 

59. 241..  4275. 163. 1.7. 954. 1.1 02. 9. 276.. 41 8. 157. 1.9. 038. 1.29, 

53.999. . 6152.97.6.1.297.5.792.4.306. .5550.97.99.5.514.5.754, 

53. 61 5..  61 97. 86. 26.. 333. 8. 689. 5. 8. .49.147.2.6.903. 1.289, 

55.6. . 49.130. . 10.. 2.844. 3. 55,. 6068. 124.7,1 .1 12,3.119, 

53. 6..  62. 105. 8.. 1692. 6. 026. 5. 4.. 53. 103 .1.3. 93 1.7. 767, 

53.97..  6459 .131. 8 ..2233. 2. 723. 3. 65.. 64. 126. 8. .6834.3.411, 

53. 118 .. 651 9. 164.9. 1.208 .1.51.2. 031. .7191. 153.1.1.362.1.950, 

5 14. 4..  3923. 152. 5 .8. 354. 2. 597. 10. 99.. 4 599. 138. 4 .4. 81 1.3. 726, 
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11 

600 

101 

10? 

103 

104 

105 

106 

111 

11? 

113 


DAT  A ( ARA (  I ) .1*781,880)/ 

$16.6, ,3773. 22*. 1,6. 28 ,.9121, 10. 5*,. 4533 ,159. 3 ,4. 832, 2. 529. 

SI  0.33, .4502, 162., 5. 132* 2.444* 10. 15 *.4471 *165. 6, 5. 378, 2. 3?R* 
$9.976, .4439, 168. *5. 72 1*2. 258, 9. 804, .4408*176.2,5.675, 1 .997, 
$14.22, .363,228.4,7.024,1.016,9.952, .4318,233.5,5.065, ,9?44, 
$9.272, .4345, 2 10., 4. 91 1, 1.258, 10.13, .4146,225.7,5.525, 1.055, 

58. 949..  4304. 2 13. 3. 5. 071. 1.221. 11. 94.. 3783. 247. 2. 6. 655,  .849, 
$8.472, .4405, 195.5,4.051 ,1.604,8.301 , .4399,203.7,3.667, 1 .459, 

$6. 567,. 4858, 193., 2. 65, 1.66, 5. 951, .5016, 196.1,2.662,1.589, 

57. 495..  4523. 251 .4 . 3 . 433 . . 86 19 . 6 . 335 . . 4825 . 255 . 1 . 2 . 8 34 . . 82?8 , 

54.3 14. . 5558.214.8.2.354.1.263.4.02. .5691 .219.9.2.402. 1.  191/ 

DAT  A  ( ARA ( I ) ,  1*881 ,1000)/ 

53. 836..  5765. 2 10. 2. 2. 742.1 .305.4.68. .5?47, 244. 7, 2. 749,  .8962, 
$3.2?3,. 5883, 232. 7, 2. 954,1 .05, 2.892,  .6204, 20  8.6,  2.415,  1.416, 

54. 728. . 5522.217. .3.091. 1.386.6.18. . 52 . 1 70 . . 4 . . 3 . 224 , 

S9.  , .47,198. .3.8,2.032,2.324, .6997,216. , 1.599,  1.399, 

51 .961 . . 7286.223. .1 .621 .1 .296. 1.75. .7427.350.1 . .9789. .5507, 

SI  0.31,. 46 13, 261. 2, 4. 7 38,. 9899, 7. 96?,.  51 9, 235. 7, 4. 34 7,  1.31 3, 

$6. 227,. 5645, 231. 9, 3. 961, 1.379, 5. 246,. 5 947, 228. 6, 4. 027,1 .43?, 

55.4  08,  .58 11, 235. 7, 3. 961, 1.358, 5. 2 18,. 5828,245., 3. 838,  1.25, 

50,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 

50,0,0,0,0,0,0,0,0,0/ 

DATA  ATM ASS/1 .008,4.003,6.939,9.012,10.811 , 12.011 , 14.007,15.999, 
A  18.998,20.183,22.989,24.312,26.982,28.086,30.974,32.064,35.453, 

R  39.948,39.102,40.08,44.956,47.90,50.942,51.996,54.9380,55.847, 

C  58.933,58.71,63.511,65.37,69.72,72.59,74.922,78.96,79.909, 83.80, 
D  85.47,87.62,88.905,91.22,92.906,95.94,99.00,101.07,102.905, 

E  106.40,107.870,112.40,114.82,118.69,121.75,127.60,126.904, 

F  131.30,132.905,137.34,138.91,140.12,140.907,144.24, 

6  145.0,150,35,151.96,157.25,158.924,162.50,164.930,167.26, 

H  168.934,173.04,174.97,178.49,180.94,183.85,186.2,190.2,192.2, 

I  195.09,196.967,200.59,204.37,207.19,208.980,242.0,210.0, 

J  222.0,223.0,227.0,227.0,232.038,231.0,238.03,237.0,242.0, 

K  243.0,245.0,249.0,249,0,254.0,252.0/ 

DATA  ERFF/0. 01 ,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1 ,0.15, 
*0.2, 0.25, 0.3, 0.35, 0.4,0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 
*0.9, 0.95,1. ,2. ,3. ,4. ,5. ,6. ,7. ,8. ,9. ,10./ 

FORMAT < I8,F14.2,3X,5(7X,A2,F10.2,  IX)  ) 

FORM AT ( 1 0 1 5 ) 

FORMAT (1H1,45X»26HENERGY  LOSS  REFERENCE  DATA// 

*30X, 10HPARTICLE  1 . 55X. 1 OHPARTICLE  2/ 

*32X,A2,1H-, I3»59X» A2» 1H-* 13// 

*5X,6HTARGET,5(5X,A2) ) 

FORMAT  <1H*,75X,5<5X.A2) ) 

FORMAT (3X, 7HMEV/ AMU, 58X,7HMEV/AMU) 

FORMAT ( 1X,F9.3,3X,5F7.3) 

FORMAT ( 1H+,65X,F9»3,3X,5F7.3) 

FORMAT (1X,F9.3,3X, 15F7.3) 

FORMAT(  1H1,45X,33HENER15Y  LOSS  VALUES  F  OR  THE  TARGET// 

*52X, 17HINC0MING  P ARTI CLE/57X , A2 , 1 H- . I 3//6X , 5HLAYER. 1517) 

FORMAT (1H1.45X.33HE ME RGY  LOSS  VALUES  FDR  THE  TARGET// 

*52X»17H OUTGOING  P ARTI CLE/57X , A2, 1 H- , 1 3//6X . 5HLA YER . 1 5 1 7 ) 

FORMAT (4X.6HENERGY) 

M  =  1 

DO  300  J*1 , 2 
DO  300  1*1,100 
DO  300  K=1 , 5 
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ACO(K»I,J)*ARA(M) 

M*M*1 

300  CONTINUE 

READ(5»600)ZIN(1) »ZIN(2) 

C  I CON*0  IMPLIES  THAT  NO  CONCEN-TABLE  NECC.  OTHERWISE  YES. 

READ (5*600) IZA, ( ZTARQ (K ) , K*1 , 5) . ICON 

C 

DO  1000  1*1 ,5 
DO  1000  J*1 » NLAYER 

CATARG( It J)*0.SATARG( I » J) *0 . SDUMT ( J ) *0 . 

1000  CONTINUE 
JZ*IZA-1 
DO  602  1*1, JZ 

READ (5,*) (ATARQ(I.K) ,K*1 .NLAYER) 

DO  602  J*1 »  NLAYER 

DUMT ( J) *DUMT( J) ♦  ATARGd,  J) 

602  CATARQ(I,J)*ATARG(I,J) 

DO  100  1*1, NLAYER 
ATARG(IZA,I)*100.-DUMT(I) 

CATARG(IZA,I)*ATARG(IZA,I) 

C  PROVISION  IS  MADE  FOR  5  CONSTITUENTS  IN  FACH  LAYER. 

ATSUM*0. 

C  ZTARG  IS  THE  ATOMIC  NUMBER  OF  EACH  CONSTITUENT. 

C  ATARG  IS  THE  ATOMIC  PERCENTAGE  OF  EACH  CONSTITUENT. 

C  THE  SUM  OF  THE  ATARGS  MUST  BE  EXACTLY  100. 

C  REL  CONCEN  IS  CALCULATED  IN  HERE. 

DO  500  lATal.5 

IF(ZTARGdAT)  .EQ.O.)GO  TO  500 
AMA*ATMASS (ZTARG d AT) ) 

I F ( AMA , FQ • 0 • ) GO  TO  500 
ATSUM*ATSUM*AMA*ATARG(IAT,I) 

500  CONTINUE 

IF  (ATSUM,EQ.O.)GO  TO  501 
DO  506  I AZ* 1 , 5 

506  CONCEN (I AZ. I) *ATARG (I AZ , I ) / ATSUM 
GO  TO  50? 

501  DO  511  I AZ*1 .5 

511  CONCEN (IAZ,I)*1. 

502  CONTINUE 

DO  505  I  AT* 1 . 5 

IF(ZTARGdAT)  ,EQ.O.)GO  TO  505 
ATARG(lAT.I)*ATARG(IAT, I) *ATMASS(ZTAR3(IAT) ) 

505  CONTINUE 

FACTOR*ATARG(  1 ,1)  ♦ ATARG  (2, 1)  ♦  ATARGd,  I)  ♦ATARG  (4, 1)  ♦  ATARG  <5.  I) 
IF  (FACTOR. LT.O.R)  GO  TO  100 
CALL  SETION(ZTARG.FIONd)  ,1) 

DO  15  K*1 ,5 

IF (ZTARG (K) .LT.l.OR.ZTARG(K) .GT.100)  ZTARG(K)=100 
15  CONTINUE 

IF ( ISSW ( 8 ) . NE. 1 ) GO  TO  3004 

WRITE (6, 11) I, FI ON (I) * ( SYMB ( ZT ARG ( K ) ) .CATARG(K.I) ,K*1,5) 

3004  CONTINUE 
100  CONTINUE 

IF (ISSW (1 ) .NE.l )GO  TO  200 

C  OUTPUT  ENFRGY  LOSS  DATA  FOR  THE  INDIVIDUAL  TARGET  CONSTITUENTS 

MASS1*IFIX (Ml*. 5) 

MASSR=IFIX (M3*. 5) 
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WRITE  (6*1 01 ) SVHB  <  Z IN  < 1 ) ) .MASS1 1 SYMB ( 2 IM < 2) ) .MASS2, ( SYMB (ZTARG ( I ) ) « 
*1*1. IZA) 

WRITE (6. 10?)  ( SYMB ( ZT  ARG ( I ) ) .1*1 .IZA) 

WRITE (6.103) 

DO  110  1*1*37 
DO  155  J*1 . 2 
IP*ZIN< J) 

PI*  ATMASS ( i ) 

IF(IP.EQ.2)PI«ATMASS(2) 

F*EREF(I)*PI/M1 
IF(J.EQ.2)E*EREF(I)*PI/H3 
DO  154  K= 1 , IZA 
JP=ZTARG(K) 

AlsACO(l.JP, IP) 

A2=AC0(2, JP.IP) 

A3=AC0(3. JP.IP) 

A4=AC0(4. JP.IP) 

A5=AC0(5. JP.IP) 

SFAC=DFAC ( JP) *.001 
SL=M*(E*1000.)  **A2 
SH=A3*ALOG( l.*A4/E*A5*E)/E 

154  SI (K,I,J)sSL*SH/(SL*SH)*SFAC 
IF (  J . EQ • 2 ) GO  TO  155 

WRITE (6.1 04 )E« (SI (K,I»1).K*1*IZA) 

155  CONTINUE 

WRITE (6. 105) E. (SI (K.I.2) .Ksl.IZA) 

110  CONTINUE 

OUTPUT  THE  ENERGY  LOSS  VALUES  FOR  THE  r I RST  15  LAYERS 
NZsNLAYER 

IF(NLAYER.GT.15)NZ=15 

DO  125  Ms  1 »  2 

IF ( M.EQ.2) GO  TO  127 

WRI TE (6 . 1 i i ) SYMB ( Z IN ( 1 ) ) .MASS1. <N.N*1.VZ) 

GO  TO  126 

127  WRITE  (6.11?)SYMB(ZIN(2) ) » MASS2 , <N.N*1.NZ) 

126  WRI TE ( 6 . 1 1 3 ) 

IP*ZIN(M) 

PI*ATMASS(1) 

IF( IP.EQ.2) PI *ATM ASS ( 2 ) 

DO  1?1  1*1,37 
FsEREF(I)*M1/PI 
IF(M.EQ.2)E*EREF( I)*M3/PI 
DO  120  Jsl.NZ 
S2( J) *0. 

DO  120  Ksl.IZA 

120  S?( J)*S2< J) *S1 (K. I.M)*CATAPQ(K, J) /100. 

WRITE(6,106)E. (S2(N) .N*1,NZ) 

121  CONTINUE 
125  CONTINUE 
200  CONTINUE 

RETURN 

FND 
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SUBROUTINE  WI DTH ( E , DE * M, I , W ID» L ) 

C  THIS  SUBROUTINE  CALCULATES  THE  STANDARD  DEVIATION  OF  THE  ENERGY 

C  LOSS  STRAGGLING. 

C  THE  EXPRESSION  FOR  ENERGY  LOSS  STRAGGLING.  WIO.  IS  FROM  THE 

C  BOHR  MODEL  AS  DESCRIBED  BY  DEARNALE Y  AND  NORTHROP.  «SEMICON- 

C  DUCTOR  COUNTERS  FOR  NUCLEAR  RADIATION*.  1963.  THE  SAME  EXPRP- 

C  SSION  IS  FOUND  IN  CHU.  MAYERtAND  N I COLET . • BACKSCATTER I NG 

C  SPECTROMETRY* ♦  1978.  THE  LINDHARD  AND  SCHARRFF • S  CORRECTION 

C  FACTOR  FOR  LOW  AND  medium  energy  projectiles  IS  USED. 

C  E  -  PROJECTILE  ENERGY  IN  MEV 

C  DE  -  ENERGY  LOSS  IN  MEV 

C  M  -  PROJECTILE  MASS 

C  I  -  IONIZATION  POTENTIAL  OF  THE  LAYER  IN  EV.  CALCULATED  IN 

C  SETION. 

C  WMAX  -  MAXIMUM  ENERGY  TRANSFERRED  TO  AN  ELECTRON  IN  EV. 

C  WID  -  STANDARD  DEVIATION  OF  THE  ENERYY  LOSS  STRAGGLING 

C  DISTRIBUTION. 

REAL  I ♦ M 

COMMON/FL AGS/FLAG. FLAGA.FLAGB.FLAGC.FLAGD.FL AGE 
COMMON/REACT/MI, M2. M3 *M4.Q,ZIN( 2) .ZTARG (5) 

COMMON/TABLES/ACO (5. 100.2) , CONCEN ( 5 .1 0 0 >  » CAT A RG ( S . 1 0 0 ) »DFAC( 100) 
COMMON/S /VALUE (2000) .ISSW(IO) 

IF  ( ISSW  (2)  .NE.DGO  TO  20 
IF(DE.LT.1.E-10)GO  TO  20 
WMAX  =  4.**E/M*5.486E*2 
Z2*0  . 

DO  1  J*  1 » 5 

1  Z2*Z2+ZTARG( J) *CATARG ( J.L) /100. 

CHI»39.5B39*E/(M*Z2) 

CLO*ALOG(WMAX/I) 

IF ( CHI-3. ) 2,2,3 

2  WID*SQRT( (*MAX*0E*l.E-6)/4.) 

return 

3  IF  (CLO.LE.O.) GO  TO  10 
WID=SQRT(WMAX/(2.*DE*1.E6*CL0)  )  *DE 
RETURN 

10  FLAGC-1. 

20  W 1 0*0  . 

RETURN 

END 
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SUBROUTINE  PICKE(F.EO.SIG.E) 


THIS  SUBROUTINE  FINDS  THE  VALUE  OF  E  SJCH  THAT  THE  INTEGRAL 
OF  THE  GAUSSIAN  DISTRIBUTION  FROM  MINUS  INFINITY 
TO  E  EQUALS  F.  THE  INTEGRAL  FOR  ALL  E  GIVES  F  =  1.00 

FAeF-0.5 
E  =  E0 

FA  IS  THE  AREA  BETWEEN  E  AND  EO. 

IF ( FA ) I0»i00»20 
F  LESS  THAN  0.5 
10  FA=0.5-F 

CALL  F INOE (FA * X) 

X=-X 

GO  TO  30 

F  GREATER  THAN  0.5 
20  CALL  FINDE(FA»X) 

30  E=E0*X«SIG 
100  CONTINUE 
RETURN 
END 


SUBROUTINE  FINDE ( F*  X ) 

X  IS  FOUND  SUCH  THAT  INTEGRAL  FROM  0  TO  X  IS  F. 

X= (E-EO) /SIG 

IF(F.LT.O.OA)  GO  TO  10 

Q  IS  THE  INTEGRAL  FROM  X  TO  INFINITY 

Q=0.5-F 

IF(O.LT.l.E-lO)  Q*1 • E-l 0 

THIS  EXPRESSION  IS  FROM  HASTINGS.  'APPROXIMATIONS  FOR  DIGITAL 
COMPUTERS'.  SHEET  88. 

U=SQRT(ALOG(l./Q**2) ) 

X=U-(2.515517+U*(0.802853+U*0. 010328) ) / 

*  (1.+U*(1.432788+U*(0.1892S9+U*0.001309) ) ) 

RETURN 

10  X=F /0 . 3R89423 
RETURN 
END 
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SUBROUTINE  DIST(EtSIG»WT t EM AX ♦ EMIN) 


THIS  ROUTINE  DISTRIBUTES  A  GAUSSIAN  INTO  'VALUE'  GIVEN  THE 
MEAN  VALUE (E) t  ANO  THE  STANDARD  DEVIATION  (SIG). 

VALUE  CONTAINS  2000  BINS. 

WT  IS  USED  TO  WEIGHT  THE  ENTIRE  DISTRIBUTION. 

EMIN  AND  EMAX  DEFINE  THE  ENERGY  SFAN  0*  VALUE. 

COMMON/S/VALUE (2000) *  ISSN (10) 

FACTOR*2000./(EMAX-EMIN) 

EA*EMIN-E 

LIMIT  THE  GAUSSIAN  TO  FIVE  STANDARD  DEVIATIONS. 

EL*E-5.*SIG 

NL*  (EL -EMIN) "FACTOR 

IF(NL.LT.l)  NL*1 

EH*F*5.*SIG 

NH*  (EH-EMIN) "FACTOR 

I F ( NH.GT, 2000 )  NH-2000 

SIGA*1./(SIG"1.4142) 

ACON*SIGA"EA 
AX«SIGA /FACTOR 
TH=  AX*  (NL-1) +ACON 
CALL  IGAUS (TH»UH) 

DO  100  I*NLtNH 

TH»AX*I>ACON 

UL«UH 

CALL  IGAUS (TH.UH) 

100  VALUE(I)«VALUE(I)>(UH-UL)*WT 
RETURN 
END 


SUBROUTINE  IGAUS(T.F) 

FIND  F8  THE  AREA  BENEATHE  THE  GAUSSIAN9  FROM  -  INFIN  TO  T 
T*(F-EO)"SQRT (0,5) /SIG 
THE  AREA  FOR  ALL  T  IS  1.0 

VcABS(T) 

Cl*0. 0705230784 
C2«0. 0422820123 
C3«  0.0092705272 
C4*  0.0001520143 
C6«  0 . 0002765672 
C5«  0.0000430638 

U«0.5/(1.*V"(C1^V*(C24'V*(C3+V*(C4^V*(C6^V*C5)  ))  )  )  )**16. 

U  IS  THE  INTEGRAL  FROM  -  INFIN  TO  -V  TAKEN  FROM  N8S-AMS-55. 

ERROR  IS  LT  l.E-7 
F«U 

IF (T.GT.0,)F»1.-U 

RETURN 

END 
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SUBROUTINE  KIN(E*ANGL»E3»XSEC»KM) 


THIS  IS  A  CLASSICAL  KINETICS  SUBROUTINE.  XSEC  IS 
THE  REACTION  IS  M1(M2*M3)M4. 

REAL  Ml  *  M2  *  M3  *  M4 

C0MM0N/REACT/M1 t M2, M3* M4.Q,ZIN( 2) »ZTARS(5) 

COMMON  /FLA0S/FLAG.FLAGA»FLAG8»FLA6C»FLAGD*FLAGE 
Ds(M1+M2)*(M3*M4) 

C=E/(D* (E*Q> ) 

A=M1»M4»C 

B=M1*M3*C 

E3s1.0*M1*0/(M2*(E+Q) ) 

C=M2*E3/D 

0=C»M4 

C=C»M3 

SET  FLAG=1  IF  DOUBLE  VALUED 
IF ( R.GT.D)  flag»i.o 
XSEC=0. 

F3  =  0. 

U=D/B-SIN(ANGL>**2 
IF (U.LT. 0 • )  RETURN 
U=SORT (U) 

E3*(E  +  Q)»B«MC0S(ANGL)*U>»*2 
V=(E3/(E*0)-B-D)/(2.*SQRT(A»C) ) 

IF(V.GT.0.7>  GO  TO  10 
IF  (V.LT.-0.7)  GO  TO  20 
ANGCMaACOS(V) 

GO  TO  30 

10  ANGCM*ASIN(SQRT(E3/(E*Q)/D)*SIN(ANGL)  ) 

GO  TO  30 

20  ANGCM*3.141592654-ASIN(SQRT(E3/(E  +  Q)/0>*SINUNGL>  ) 
GET  CENTER  OF  MASS  CROSS-SECTION. 

0  CALL  SIG(E»ANGCM,XSEC.KM) 

CONVERT  TO  LABORATORY  CROSS-SECTION. 
XSEC=XSEC*E3/(SQRT (A*C)*U*(E*Q)  ) 

RETURN 

END 


LABORATORY. 


45 


no  o  noon 


SUBROUTINE  SIQ(E* AN6CM*XSEC*KM) 

THIS  SUBROUTINE  CALCULATES  THE  REACTION  XSEC  FOR 
M1(M2»M3)M4  RUTHERFORD  IF  M2*56  »  FOR  0(0»P1)  IF  M2«16»  AND  FOR 
N(D,P5)  FOR  M2*14. 

FOR  C(D.PO)  FOR  M2*12. 

REAL  M1,M2,M3»M4 

COMMON/REACT/MI ,M2,M3, H4»Q»ZIN (2) ,  ZTARG(5) 

DIMENSION  EN(40) *XO(40) *  XN ( 40 ) »XC(40) 

INTEGER  ZIN*ZT AMG 

DATA  (EN(I) .1*1 .40) /.3, .31t.32t.33. ,34. .35, .36, .37,. 38, .39* .40, 
$»41»,42»*43».44».45».46»»47»*48»,49».50».51».52».53»,54..55,.b6, 
$.57*.58»,59».60».61».62».63».64».65*.6S».67*.68».69/ 

DAT  A  ( XO  ( I )  *I*1»40)/0.».01,.025*.03*.035».04,.046*.055*.066,.082, 
$.10, ,129* ,160*. 200, .235*. 270, .310, *350, .390, .440, .500, .580, .690* 
$.80  0,.  930, 1.050, 1.21 0*1. 380, 1.51 0*1. 64 0,1. 720, 1.785, 1.835, 1.875, 
$1.900*1.920*1. 940*1. 955*1.965*1. 970/ 

D AT  A ( XC ( I ) ,  1*1.40 )/0. 208, 0.20  9* 0.21 0*0. 21 1*0. 21 5. 0.21 9, 0.220,0.??, 
$0.224,0.225*0.229,0.230,0.231,0.240,0.249,0.251,0.261,0.271  , 
$0.288,0.305,0.324,0.340,0.360*0.385,0.411,0.445,0.480,0.520,0.560, 
$0,610*0.660*0.715,0.772,0.840,0.915,0.999,1.100,1.250,1.430,1.610/ 
D AT  A ( XN ( I ) , I *1,40)/. 01 5,. 020*. 024*. 028,. 037,. 047,. 058,. 074,. 092,. 1 
*14,. ]39,. 169*. 202 *.235*. 279,. 31 8,. 368 *.423, .4 78,. 539, . 599 , . 66 0 , . 72 
*7,  .793,  .849, .905, .939, .961  , . 95 1 , . 929 , . 897 , . 886 , . 887 , . 898 , . 954 , 1 . 0? 
*7,1.111,1.190*1.268,1.369/ 

Z1*ZIN(1) 

Z2*ZTARG (KM) 

IF  Q.EO.O.  ELASTIC  SCATTERING 
IF(Q.EQ.0.)GO  TO  10 
CM2=M2-16. 

CM2.ABS (CM?) 

I  CMs  1 

IF(CM2,LT. 1.0)60  TO  20 
CM3.M2-14. 

CM3.ABS (CM3) 

I  CM*2 

I F ( CM3.LT .1.0)60  TO  20 
CM4*M2-12. 

CM4*ABS (CM4) 

I  CM*3 

I F ( C84 . LT . 1 . 0 ) GO  TO  20 

THIS  SECTION  CALCULATES  THE  RUTHERFORD  CROSS  SECTION  OF  THE 
REACTION  M2(M1,M3)M4 
10  CS=]./SIN(ANGCM/2.) 

CS4=CS**4 

XSEC* 1  .296* (Z1*Z2/E) **2  *( ( Ml ♦ M2  I /82 ) **2*CS4 
RETURN 

20  IF (E.LT.EN(l) . OR . E . GT . EN ( 40 ) ) GO  TO  100 

DO  21  1*1,39 

EHI=EN(IM) 

ELO=EN( I) 

IF(E.GE.ELO.AND.E.LE.EHI)GO  TO  22 

21  CONTINUE 
?2  CONTINUE 

DEE  *EHI -ELO 
DE*E-ELO 

IF(ICM-?)?3,24,25 
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THIS  SECTION  CALCULATES  THE  0(0»P1)  CROSS  SECTION  FROM 
DATA  TABLE  FI 39  CROSS  SECTION  VALUES  ARE  FROM  LONGEQUE  ET . AL • * 
JOURNAL  DE  PHYSIQUE84UNE  65.  THEY  ARE  GOOD  ONLY  FOR  160  DEGR  LAR 
ANGLE 

XSEC*X0(I) ♦ <  XO ( I ♦ 1 ) -XO ( I ) ) *DE/DEE 
RETURN 

CROSS  SECTION, VALUES  ARE  EXTRAPOLATED  FROM  THE  OATA  OF  NIILER  ANn 
BIRKMIRE (REF  TO  BOOTH  ET.AL.  (PROC  PHYS  SOC. 71 1 ( 1 957 ) , 325 . ) 
XSEC«XN(I)*(XN(l4l)-XN(I) )*DE/0EE 
RETURN 

THIS  SECTION  CALCULATES  THE  C(D.PO)  CROSS  SECTION  AND  IT  IS  GOOn 
FOR  160  DEGR  LAB  ANGLE 
C(D.PO)  HUEZ  ETAL.  NIM  105(1972)197 
25  XSEC*XC(I) ♦ (XC(I*1)”XC(I) )*DE/DEE 

RETURN 

C  ENERGY  OUT  OF  RANGE  OF  TABLE  INPUT 

100  XSEC=0. 

RETURN 

C  MASS  OF  TARGET  IS  rfRONGBRESTART  PROGRAM 

END 
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SUBROUTINE  SHAPE  ( E . DE t W I D. FMASS* ION. FI  ON , THI CK . THK VAR ♦ LA YER ) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  MEAN  ENERGY  AND  STANDARD  DEVIATION 

C  AS  PARTICLES  TRAVERSE  ALL  BUT  THE  FINAL  LAYER. 

C 

C  E  IS  THE  INITIAL  PARTICLE  ENERGY. 

C  FMASS  IS  THE  PARTICLE  MASS. 

C  THICK  IS  THE  ARRAY  OF  LAYER  THICKNESSES. 

C  F ION  IS  THE  ARRAY  OF  IONIZATION  POTENTIALS. 

C  LAYER  IS  THE  NUMBER  OF  THE  DESIRED  LAYER. 

C  ION  IS  THE  ID  FOR  THE  PARTICLE*  1  IS  BEAM  AND  2  IS  OUTGOING. 

C 

COMMON/SHA/NSLA(100) • MT ( 1 00 ) 

DIMENSION  THICK ( 1) ♦ FI ON(l) . THKVAR ( 1 ) 

C  CUTOFF  SET  AT  1  KEV 

IF(E.LT.l.E-3)  GO  TO  20 
WID  =  0. 

DE=WID2=0.0 
IF (LAYER, LE, 1 )  RETURN 
EA  =  F 

NLAYcLAYER-1 

IF ( LAYER. GT.2.AND.THKVAR(15 .NE. 0 . AND . I DN. EQ . 2 >  NL A Y*L A YER-2 
00  10  I « 1 ♦ NL AY 
1 1  =  T 

IF ( ION.EQ,?)  II»LAYER-I 
IF (THICK ( II ) .LE.O.O)  GO  TO  10 
UA  =  1. 

N5L AR*NSL A (II) 

THKN«THICK(II)/NSLAB 
DO  30  J*1  * NSLAB 
31  THKAeTHKN*UA 

CALL  ELOSS ( EA* 1 1  *  ION*  THKA* OE) 

IF (EA.LT ,DF)  GO  TO  20 

CALL  WIDTH (EA*DE»FMASS.FION (II ) »WID»LAYER) 

WID?=WID2+WID**2 

ea=ea-de 

30  CONTINUE 
DE=F-EA 

10  CONTINUE 

W I D*FQRT ( WI 02) 

RETURN 

FITMER  OE  CANNOT  BE  FOUND  OP  PARTICLE  IS  STOPPED 
IN  EITHER  CASE*  PARTICLE  IS  ASSUMED  STOP»ED 
DE  GREATER  THAN  E  IS  USED  TO  FLAG  THE  CALLING  PROGRAM. 

0  0  E  =  E  ♦  E 

WID*0.0 
RETURN 
END 
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SUBROUTINE  ELOSS ( E , LAY , I  ON , T , OE ) 

c  this  routine  calculates  the  energy  loss  of  protons,  deuterons, 

C  TRITONS  ,  ALPHAS  OR  HE3»S  AT  THE  ENERGY  E  (MEV).  Ml  AND  ZlN(l) 

C  ARE  THE  MASS  AND  CHARGE  OE  THE  INCIDENT  PARTICLE,  M3  AND  ZIN(?) 

C  ARE  THE  MASS  AND  CHARGE  OE  THE  OUTGOING  PARTICLE.  ZTARG(I) 

C  GIVES  THE  CHARGE  OE  THE  I'TH  TARGET  CONSTITUENT,  CATARG(I,J) 

C  GIVES  THE  ATOMIC  PERCENT  OE  THE  I*TH  CONSTITUENT  IN  THE  J*TH 

C  LAYER.  T  IS  THE  SLAB  THICKNESS.  THE  EXPRESSION  EOR  DE  IS  FROM 

C  A/Z.  FOR  ALPHAS,  IT  IS  GOOD  10  TO  10,000  KEV  AND  FOR  PROTONS, 

C  GOOD  FROM  10  TO  1000  KEV.  EOR  PROTONS,  IE  USEO  AT  10  MEV,  Thf 

C  ERROR  IS  LESS  THAN  5  PERCENT. 

COMMON/TABLES/ACO (5,100,2) , CONCEN ( 5 , 1 0 0 ) ,CATARG(5,100),DFAC(100) 
COMMON/REACT/MI ,M2,M3,M4,Q,ZIN(2) »ZTARS(5) 

common/elo/iza 

REAL  M1,M?,M3,M4 
integer  ZIN,ZTARG 
IF (E.EQ.O. ) GO  TO  11 
IP=ZIN( ION) 

OE=  0  . 

PI=1 .00B 

IF  (IP.EQ.?)PI*4.003 
EC=E*PI/M1 

IF  (I0N,EQ.?)EC=E*°I/M3 
DO  10  J=1,IZA 
JP  =  ZT  ARG ( J ) 

IF  (  JP.EO.O) GO  TO  10 
A1=AC0(1,JP,IP) 

A2=AC0(2, JP,IP) 

A3=AC0(3, JP,IP) 

A4=AC0 (4, JP, IP) 

A5=AC0(5, JP,IP) 

SF AC=DF  AC ( JP ) *0.001 

SL=A1*(EC*1000.)**A2 

SH=A3*AL0G(1,*A4/EC*A5*EC)/EC 

OE  =  DE*CATARG(J,LAY)*SL*SH*SFAC/(100.*(SL*SH)  ) 

10  CONTINUE 

de=de*t 

11  EL=.01 
EH=1 0 • 0 

IF (EC.LT.EL.OR.EC.GT.EH) WRITE (6, 1 ) EC 
1  FORMAT  (4X,6>HENERGY,F10. 4, 46HMEV.  OUT  OF  RANGE  OF  ANALYTICAL  A/Z  fx 

*PRESSION) 

RETURN 

END 
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SUBROUTINE  OUTPUT ( YMAX2 ) 


THIS  ROUTINE  PREPARES  THE  CALCULATED  AND  EXPERIMENTAL  DISTRI¬ 
BUTIONS  FOR  PLOTTING.  IT  DOES  LEAST  SQUAREO  FITTING  OF  THE  TWO. 

COMMON/OU/XUIO)  »Y(410),YY(410)»XX(410) 

COMMON/S/VALUE (2000) tlSSW(iO) 

COMMON/OEX/EMIN.EMAX 

C0MM0N/REACT/M1,M2,M3*M4,Q,ZIN(2) ,ZTARG(5) 

REAL  Ml »  M2  *  M3 » M4 
COMMON/PLO/LP(400»125) 

DIMENSION  AX(400)*AY(400),AYX(400)»EXX(400) 

CONDENSE  THE  OUTPUT  BY  A  FACTOR  OF  FIVE. 

00  1111  1.1,410 
X(I).0.SY(I)«0.SYY(I)*0. 

Ill  CONTINUE 

DELE. (EMAX-EMIN)*5./2000. 

YMAXrO. 0 

READ(5,111)XSC,YSFA 
READ (5*111) AY1,AY2»BY1»BY2 
DO  10  1*1 ,400 
X ( I ) «I 
J=I»5-4 

Y  ( I )  .VALUE  (J)  ♦  VALUE  (J*l  )  +VALUE  ( J*2)  ♦  V  ALUE  ( J*3 )  ♦  VALUE  (  J* 4 ) 

111  FORMAT ( 4F 1 0 • 5 ) 

10  CONTINUE 

IF ( ISSW (7) .EO.l) WRITE (6,75) 

DO  151  1*1,400 
XX(I).X(I) 

151  YY(I).Yd) 

15  CONTINUF 

ZERO  THE  ARRAY  USED  FOR  ACCUMULATING  THE  DISTRIBUTION. 

ENTRY  OUT  1 
DO  100  1=1,2000 
100  V  ALUE ( I ) =0 , 0 
RETURN 

PLOT  THE  DATA  FOR  COMPARISON. 

ENTRY  OUT? 

READ (5, 1000) BIAS, CALX 
000  FORMAT (FB,0»1PE12.5) 

BIAS  IS  THE  NO  OF  CHANNELS  THE  DATA  ARE  OFFSET  FROM  ZERO 
CALX  IS  THE  ENERGY  SCALE  CALIBRATION  IN  MEV/CHANNEL 
1  FORMAT (10FB.0) 

READ(5,11?) ICNT 

C  ICNT  GIVES  THE  NUMBER  OF  DATA  POINTS  TO  RE  READ.  ICNT. LE. 400 

READ (5, 21) ( X ( I ) ,Y(I) ,1*1, ICNT) 

SUMXC»SUMX?*0. 

YMAX.O. 0 
NPTS=0 

YMAX?*0SIX«0SKSS*0 
DO  ?0  1=1 , ICNT 
EX= ( X (I) -BI AS)*CALX 
DO  23  JC=1 ,400 
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EC=EMIN*DELE*XX ( JC) 

IF(EC.GT.EX.AND.EX,GE.EMIN)GO  TO  24 

23  CONTINUE 

24  AY(I)*(YY( JC)-YY( JC-1) ) * ( EC-EX ) /DELE* Y Y ( JC- 1 ) 
IF(EX.GE.EMIN.AND.IX.EQ.O) ix«i 
IF(EX.GE.AY1.AND.EX.LE.AY2)G0  TO  999 

GO  TO  998 

999  SUMXC*SUMXC*AY ( I ) *Y ( I ) 

SUMX2*SUMX?*Y ( I )*Y( I) 

99fl  CONTINUE 

AX ( I ) *X ( I ) 

AYX ( I ) *Y  ( I ) 

IF(EX.LE.EMAX)NPTS«NPTS*1 
20  CONTINUE 
25  CONTINUE 

ASUM=SUMXC/SUMX2 
NPSsNPTS 
DO  SO  J*1 *  I CNT 
Y(J)=AYX(J)*ASUM 
50  YMAXsAMAXl (YMAX»Y( J) *AY( J) ) 

YSC=YMAX/YSF  A 

CALCULATION  OF  THE  XISQ  OVER  ENERGY  REGION  AYI  TO  A Y2 

X  I SQ  =  0  # 

K  S=  0 

5  FORMAT(lHl) 

NK  =  0 

DO  980  J*1 » I CNT 

AYX?r0  %  AY22*0 
EX  =  (  AX (J) -BIAS) *CALX 
EX  1 =FX 

CARFFUL***  THE  FOLLOWING  STATEMENT  WITH  ITS  NUMBERS  IS  GOOD  ONLY 
FOR  C » N » 0  ANALYSIS  AT  670  KEV. 

IF(EX.GT.?.5)EXl«EX-0.75 
AYX( J) =ASUM*AYX ( J) /YSC 
AY ( J) *AY ( J) /Y  SC 

IF(NK.EQ.J) AYX2=( AYX ( J-l) *AYX ( J) ) /2. 

IF(NK.EQ.I) AY22*(AY( J-1)*AY(J) )/2. 

NK=NK*1 

IF ( NK ,EQ#  2) NK*0 

IF ( ISSW (7) . NE. 1 ) GO  TO  2000 

IF(EX.GT.EMAX)GO  TO  2000 

WRITE (6» 1 12) KS,EX. AYX (J) .AY (J) ♦ EX1» AYX2.AY2? 

112  FORMAT ( 15.2 ( 1  OX. F3.4.2F10. 3) ) 

2000  CONTINUE 
KSS=KSS+1 
EXX ( KSS) =EX 
YEXsAYX ( J) 

IFIYEX.EQ.O.) yex*i. 

996  IF(FX.GE.AY1.AND.EX.LE.AY2)G0  TO  991 
GO  TO  980 

991  XISQ=XISQ+YSC*(AY (J)-AYX(J) ) **2/ ( YEX» ASUM ) 

KS=KS*1 

980  CONTINUE 

SK  =  KS  5  X I SQ=X I SQ/SK 
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C  PLOT  THE  EXPERIMENTAL  DATA 

WRITE(6t2)CALXfBIAS»YMAX»XISQ,SKf ASUM 
2  F0RM*T(//1AH  PLOT  EXP/CAL  0ATA/t2X»25HTHE  ENERGY  CALIBRATION  IS, 

*F10.6»36H  KEV/CH  AND  THE  ZERO  LEVEL  BIAS  IS,F10.2,9H  CHANNFLS/ 
*21H  FULL  SCALE  VALUE  IS  » 1PE12.5,/20H  XISQ  PER  POINT  IS  tlPE12.St 
*6H  WITH  » 0PF10,2,15HPOINTS  INCLUOEO/12H  ASUH  IS  ,1PE12.5> 
IF(ISSW(5) ,NE.1)00  TO  1515 
CALL  PLOT1 (AYX,AY»EXXfYSC»NPSf IXtKSSflCNT) 

ISIS  CONTINUE 
RETURN 
END 
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SUBROUTINE  ®LOT 1 ( AYX , AY , EEX . YSC » NPS , I  X , KS , IC NT ) 
COMMON/PLO/LPUOO»125) 

DIMENSION  AYX ( 1 ) «  A Y ( 1 ) • EEX ( 1 ) 

DIMENSION  I S ( 4) 

DAT  A ( IS ( | ) *  1*1  *  4) /1H  .lHO.lH.tlHS/ 

NPTS*NPS-1 
NPTS*IX»NPTS 
DO  1  1*1,125 
DO  1  J*1 . ICNT 
1  LP(J.I)*IS(1) 

IXX=IX 

KL*NPS 

C 

C  THIS  PLOTS  THE  CALCULATED  DATA 

C  «•«•«•«•*«•******###***##*#***#*#***#######» 

DO  5  K«  IX.NPTS 
XIsAY  < K ) 

IXI«XISXII*IXISDX*XI-XIISIF(DX.GE.0.5) I X  I  =  I X  I ♦ 1  % J*  I  X  I 

IF(J.LE.O) J*1 

LP(KL.J)*IS(3) 

C 

C  THIS  PLOTS  THE  EXPERIMENTAL  DATA 

C 

X 1  =  AYX ( K) SIXI*XI$XII*IXISDX*XI-XIISIF OX.GE.0.5) I X I  *1 X I  ♦  1  *  I  =1  X I 
IF(I.LE. 0)1*1 
LP ( KL* I ) *IS (2) 

C 

c  a*********#***#.*****#****************** 

IF(J.EQ.I)LP(KL,I)*IS(4) 

KL*KL-1 

IF (KL.LE.O) KL*1 
5  CONTINUE 

KM*1SKL*NPSSN*9 
WRITE (6,29) 

?9  FORMAT  <  1 H 1 ) 

WRITE  (6,31) 

31  FORMAT ( 1H  .50X, 13HPR0FILE  CURVE) 

DO  50  K*1»NPS 

IF (KM.EQ.K. AND. EEX (KM) .NE.O) GO  TO  40 
WRITF  (6.30)  (LP(KL.I )» 1  =  1.120) 

30  FORMAT(7X.1HI,120A1.1HI) 

KL*KL-1 

IF (KL.EQ.O) KL*1 
GO  TO  50 

40  WRITE (6,35) EEX (KM) , (LP(KL.I) ,1*1,120) 

35  FORMAT (1X.F5.2, IX »1HI.120A1,1HI) 

KM*KM*N 
N=  1  0 

KL*KL-1 

IF(KL.EQ.0)KL=1 
*0  CONTINUE 

WRITE (6.51) 

51  FORMAT(8X,12(10H  + - ),2H*0 

WRITE (6,55) 

55  FORmaTUH  ,  50  X ,  5HY I  ELD ) 

IX=IXX 
RETURN 
END 
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SUBROUTINE  SWITCH 

COM  MON/ S/V  ALU  E  ( 2000  )  .ISSWUO) 

C  A  SWITCH  IS  DESIGNATED  BY  ISSW(I)  IN  THE  PROGRAM. 

C  THE  SWITCHES  ARE  BEING  USED  AS  FOLLOWS. 

f  1  -  OUTPUTS  ENERGY  LOSS  INFORMATION.  BOTH  FOR  TARGET  CONSTITUENTS 

C  (ELEMENTS)  AND  COMPOSITE  LAYERS. 

C  2  -  WHEN  ON  INCLUDES  ENERGY  LOSS  STRAGGLING  IN  THE  CALCULATION 

C  3  -  OUTPUTS  THE  PR  INFO,  USEFUL  THINGS  ARE  INCIDENT.  SCATTERED 

C  AND  OUTGOING  ENERGY  AND  CROSS  SECTIONS. 

C  CURRENTLY  NO-OP 

C  4  -  OUTPUTS  ALL  PERTINENT  INFORMATION  ON  THE  RUN  BEING  DONE  - 

C  REACTION.  GEOMETRIES.  ETC,  SHOULD  ALWAYS  BE  ON, 

C  5  -  GIVES  THE  PRONTER  PLOT  OF  EXPERIMENTAL  AND  CALC.  SPECTRA 

C  6  -  OUTPUTS  ENERGIES.  CROSS  SECTIONS  AS  A  FUNCTION  OF  SLAB. 

C  VERY  USEFUL  WHEN  ITERATING  CONCENTRATIONS  TO  FIT  DATA.  CAN 

C  BE  A  VERY  LONG  OUTPUT. 

C  7  -  OUTPUTS  NUMERICAL  LIST  OF  CHANNEL  NO.  VS  EXPERIMENTAL  AND 
C  CALCULATED  YIELDS. 

C  B  -  OUTPUTS  TARGET  COMPOSITION  DATA. 

C  9  -  PRINTS  THIS  SWITCH  INFORMATION. 

C  10  -  OUTPUTS  CONCENTRATION  TABLE. 

READ  (5.2) (ISSW(I) .1-1.10) 

?  FORMAT ( 1 0 1 1 ) 

IF ( I SSW (9) ,NE.1)G0  TO  5 
WRITE (6. 10) 

10  FORMAT ( l H 1 1 

WRITE (6. 11 ) 

11  F0RMAT(1X,»  switch  1  ON  FOR  ENERGY  LOSS  INFORMATION*/ 

*1X.»  SWITCH  2  IS  ON  FOR  ENERGY  LOSS  STRAGGLING  CALCULATION*/ 

*1X.*  SWITCH  3  ON  FOR  PUBLIC  RELATIONS  INFORMATION*/ 

*1X.*  SWITCH  4  ON  FOR  PERTINENT  DATA  ON  CURRENT  CALCULATION*/ 

*1X.*  SWITCH  5  ON  FOR  PRINTER  PLOT  OF  RESULTS*/ 

*1X.*  SWITCH  ft  ON  FOR  USEFUL  DATA  WHILE  DOING  I NTERAT I ONS * / 

*1X.*  SWITCH  7  ON  FOR  NUMERICAL  LIST  0*  PLOTTED  RESULTS*/ 

*1X.*  SWITCH  8  ON  FOR  TARGET  COMPOSITION  DATA*/ 

*1X,*  SWITCH  9  ON  FOR  SWITCH  INFORMATION  TABLE*/ 

*1X,*  SWITCH  10  ON  FOR  LAYER  CONCENTRATION  TABLE  * / 

*1X.*  NORMAL  RUNNING  SHOULD  HAVE  SWITCHES  4.5.7.8.10  ON*) 

5  RETURN 
END 
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o  o  o  n  o  o 


SUBROUTINE  CONCTAB ( NLAYER , THI CK ) 

COMMON/TABLES/ACO (5*100*2) * CONCEN ( 5 « 1 00 ) ,CaTaRG(5.100),DFAC(100) 
COMMON/TSYMB/SYMB( 100) ♦ ICON 
COM MON/RE ACT/M 1 »M2»M3»M4»Q»ZIN(2) ,ZTAR3(5) 

DIMENSION  Cl (5) ,CJ(5> 

INTFGER  zin.ztarg.symb 
DIMENSION  THICK ( 1 ) 

THIS  BLOCK  OF  WRITE  STATEMENTS  ARE  USED  ONLY  IF  THE  CONCENTRATIONS 
OF  THE  ELEMENTS  ARE  REQUIRED.  CONCENTRATIONS  ARE  BEING  CALCULaTfh 
IN  THE  MAIN  PROGRAM.  SEE  STATEMENT  600. 

THIS  ROUTINE  CALCULATES  THE  CONSTANT  FOR  EACH  LAYER  FOR  EACH  ElT 
THE  FORMULA  ISCONSTANT  *  CONCEN  *  A VO 3R ADO  NUMBER  *  THICKNESS 
AVOGRADO  NUMBER  EQUAL  6.025*10**23. 

I F ( I  CON . EQ. 0 ) GO  TO  4005 
WRITE (6. 1000) 

1000  FORMAT ( 1 H 1 ) 

WRITF (6.4000) 

40  0  0  FORMAT ( ?0 X. 10  OH********************** **********************  ******* 
******************************#*******************) 

WRITE (6.2000) 

WRITE (6.4001) 

4001  FORMAT (60X.40HCONCENTRATION  OF  ELEMENTS ( AT/CM2) /LAYER  ) 

WRITE (6,2000) 

WRITE  (6.4000) 

WRITE (6.2000) 

WRITE (6,4002)  ( SYMrt ( ZT  ARG ( K ) ) ,K  =  1.5) 

4002  F ORMAT ( 20H  LAYER  THI CK ( MG/CM2 ) , 7X , A2 , ?2X , A2 , 2 1 X , A2 . 21 X , A 2 . 24X , A? ) 
WRITE(6.2000) 

WRITE(6.4003) 

4  003  FORMAT (20X,1X,5HCI(1) ,7X,5HCJ(1)  .7X.5HCI (?) ,7X,SHCJ(2) ,7X,5HCI(3)  . 
*7X  »  SHC J ( 3 ) .7X.5HCI (4)  »7X.5HCJ(4) .7X.5HCI (5) ,7X,5HCJ(S) ) 

WRITE  (6,2000) 

2000  FORMAT ( /) 

DO  1  1*1,5 
Cl (  I)*0. 

1  CJ(I)*0. 

A V0  =  6 » 025  *  (10. **23) 

DO  3  J*l. NLAYER 
DO  2  1*1.5 

Cl ( I)*CONCFN( I.J) *AVO*THICK ( J)*.001 

2  CJ ( I ) *CJ ( I ) *CI ( I ) 

WRITE (6.4)  J, THICK ( J) » (Cl (K) ,CJ(K) *K*1.5) 

4  FORMAT (1X,I4,3X,F4.2,4X, (5(2(1PE10.3«2X) ) ) ) 

3  CONTINUE 
4005  CONTINUE 

RETURN 

END 
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APPENDIX  B  -  PROGRAM  INPUTS 


All  inputs  required  to  run  PROFILE  are  listed  in  this  appendix. 
Cards  marked  with  *  may  require  more  than  a  single  card. 


CARD 

1 

2 


3 


4* 


5 


6 


7 


FIELD 

(20A4) 

(1011) 


(8F10.0) 


(8F10.0) 

(8510.0) 


(8F10.0) 


(8F10.0) 


DESCRIPTION 


TITLE.  A  descriptive  title  for  the 
calculation  being  done. 

ISW(I) ,1=1, 10.  Sense  Switches  used  to 
control  calculations  and  outputs.  0  is 
off,  1  is  on.  See  the  listing  of  Sub¬ 
routine  Switch  for  description  of  each 
switch  function. 

CHAR,  OMEGA.  CHAR  is  the  integrated 
beam  change  and  OMEGA  is  the  detector 
solid  angle.  These  values  are  used  if 
absolute  fits  are  desired. 

Ml , M2 , M3 , M4 , QREACT , EXC I T 

Blank . 

This  set  of  cards  is  used  to  define  the 
reaction  between  the  beam  and  target. 
There  is  a  different  4  card  for  each 
target  constituent,  up  to  a  maximum  of 
five.  A  blank  card  ends  the  list  of 
reaction  cards.  The  reaction  is  defined 
as  M2 (Ml, M3) M4  with  QREACT  being  the 
reaction  Q-value  and  EXCIT  being  the 
energy  of  excitation  of  a  state  in  M4. 

E , THETAT , THETAD , EMAX , EMI N , OTHET .  This 
card  defines  the  beam  energy  and  target 
and  detector  geometries.  E  is  the  beam 
energy,  EMAX  and  EMIN  are  the  energy 
limits  in  MeV,  over  which  the  calcu¬ 
lation  is  done,  THETAD  is  the  detector 
angle  between  target  normal  and  beam 
and  DTHET  allows  for  a  non-flat  target 
surface  where  the  beam  hits  it. 

RS,DS,RT,RD  and  DD. 

RS  and  DS  define  beam  divergence,  RD 
and  DD  define  detector  opening  and 
RT  is  beam  spot  radius  at  the  target. 
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CARD  FIELD 

DESCRIPTION 

8  (1015) 

NGS,NGT,NGD 

Number  of  points  used  at  the  slit  (NGS) , 
target  (NGT)  and  detector  (NGD)  to  do 
the  calculation.  All  are  set  to  1  unless 
large  beam  spots  or  detector  apertures 
are  used. 

9  (8F10.0) 

DETRES,ERES. 

Resolution  of  the  detector  (DETRES)  and 
beam  energy  (ERES)  in  keV. 

10  (1015) 

NLAYER, NSLAB , NDI V . 

NLAYER  is  the  number  of  layers  into 
which  the  target  is  divided,  must  be 
<  100.  NSLAB  is  the  number  of  slabs  into 
which  the  first  layer  is  further  subdivided 
All  layers  are  subdivided  into  the  same 
thickness  as  the  slabs  of  the  first  layer. 
NDIV  is  the  number  of  points  used  to 
represent  the  outgoing  particle  Gaussian. 

Inputs  11  through  13  are  in  Subroutine  SETTAB. 


11  (1015) 

ZIN(l) ,  ZIN (2) 

ZIN’s  are  the  atomic  numbers  of  the 
incident  (1)  and  outgoing  (2)  particles. 

12  (1015) 

IZA,  ( ZTARG (I), 1=1, 5) 

IZA  is  the  number  of  target  constituents 

The  ZTARG* s  are  the  atomic  numbers  of 
the  constituents. 

13*  Unformatted 

ATARG(I,K),K=1, NLAYER 

ATARG(I,K)  is  the  relative  percentage 
of  the  element  I  in  layer  K.  A  new 
card  must  be  used  to  begin  the  list 
for  each  element.  I  <  5. 

14*  Unformatted 

THICK  (I),  1=1,  NLAYER 

THICK* s  are  the  thickness  (in  mg/cm  ) 
of  the  layers. 

Inputs  15  through  19  are  in  subroutine  OUTPUT 


15  (4F10.5) 

XSC, YSFA 

XSC  is  currently  non-operative 

YSFA  is  the  Y-scale  factor  for  printer 
plot,  must  be  <  120. 
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CARD 

FIELD 

DESCRIPTION 

16 

(4F10.5) 

AY  1 ,  AY  2 

AY1  and  AY2  are  the  lower  and  upper 
energy  limits  over  which  the  least 
squares  fitting  is  done. 

17 

(F8 .0, 1PER. 5) 

BIAS, CALX 

Experimental  energy  spectrum  parameters 
such  that 

E=*X-Bias) *CALX,  where  E  is  the  energy 
and  X  is  the  experimental  channel 
number. 

18 

(15) 

ICNT 

ICNT  gives  the  number  of  experimental 
data  points  to  be  input.  ICNT  <  400. 

19* 

(10F8.0) 

(X(I),Y(I), 1=1, ICNT) 

X  and  Y  are  the  channel  number  and 
number  of  counts  in  that  channel, 
respectively. 
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